Removing useless files
authorPeter Hristov <phristov@pb-d-128-141-110-112.cern.ch>
Thu, 5 Dec 2013 14:52:26 +0000 (15:52 +0100)
committerPeter Hristov <phristov@pb-d-128-141-110-112.cern.ch>
Thu, 5 Dec 2013 14:52:26 +0000 (15:52 +0100)
22 files changed:
STARLIGHT/starlight/TStarLight/TStarLight.h~ [deleted file]
STARLIGHT/starlight/TStarLight/testsl.C~ [deleted file]
STARLIGHT/starlight/include/inputParameters.h~ [deleted file]
STARLIGHT/starlight/include/inputParser.h~ [deleted file]
STARLIGHT/starlight/include/nBodyPhaseSpaceGen.h~ [deleted file]
STARLIGHT/starlight/include/singleton.h~ [deleted file]
STARLIGHT/starlight/include/slmutex.h~ [deleted file]
STARLIGHT/starlight/include/starlightStandalone.h~ [deleted file]
STARLIGHT/starlight/src/beam.cpp~ [deleted file]
STARLIGHT/starlight/src/beambeamsystem.cpp~ [deleted file]
STARLIGHT/starlight/src/gammagammaleptonpair.cpp~ [deleted file]
STARLIGHT/starlight/src/inputParameters.cpp~ [deleted file]
STARLIGHT/starlight/src/inputParser.cpp~ [deleted file]
STARLIGHT/starlight/src/nBodyPhaseSpaceGen.cpp~ [deleted file]
STARLIGHT/starlight/src/nucleus.cpp~ [deleted file]
STARLIGHT/starlight/src/wideResonanceCrossSection.cpp~ [deleted file]
STARLIGHT/test/Config.C~ [deleted file]
STARLIGHT/test/Config_.C~ [deleted file]
STARLIGHT/test/galice.root [deleted file]
STARLIGHT/test/rec.C~ [deleted file]
STARLIGHT/test/sim.C~ [deleted file]
STARLIGHT/test/syswatch.log [deleted file]

diff --git a/STARLIGHT/starlight/TStarLight/TStarLight.h~ b/STARLIGHT/starlight/TStarLight/TStarLight.h~
deleted file mode 100644 (file)
index 66731dd..0000000
+++ /dev/null
@@ -1,131 +0,0 @@
-////////////////////////////////////////////////////////////////////////
-//
-// Copyright 2013
-//
-////////////////////////////////////////////////////////////////////////
-//
-// File and Version Information:
-// $Rev::                         $: revision of last commit
-// $Authro::                      $: Author of last commit
-// $Date::                        $: Date of last commit
-//
-// Description:
-//     TStarLight.h is the include file defining data members and
-// functions specifically needed to implement an interface of STARlight
-// to ROOT's TGenerator (a part of ROOT's Virtual Monte Carlo).
-//
-// Based on work done by Bjoern Nielson
-////////////////////////////////////////////////////////////////////////
-#ifndef STARLIGHT_TSTARLIGHT_H
-#define STARLIGHT_TSTARLIGHT_H
-#include <cstring>
-
-#include <TGenerator.h>
-#include <TString.h>
-#include <TMath.h>
-
-#ifdef __CINT__
-#  undef __GNUC__
-#  define _SYS__SELECT_H_
-struct pthread_cond_t;
-struct pthread_mutex_t;
-#endif
-
-// Star Light includes
-#include "starlight.h"
-#include "upcevent.h"
-#include "inputParameters.h"
-#include "randomgenerator.h"
-
-// 
-class TObjArray;
-class TClonesArray;
-
-class TStarLight : public TGenerator {
- public:
-  TStarLight();
-
-  TStarLight(const char* name,
-            const char* title,
-            const char* slConfigFile="./slight.in");
-
-  virtual ~TStarLight();
-
-  // TGenerator overloaded methods
-  virtual void       GenerateEvent();
-  virtual Int_t      ImportParticles(TClonesArray *particles, Option_t *opt="");
-  virtual TObjArray* ImportParticles(Option_t *opt="");
-  virtual void       SetParameter(const char* key, Double_t val);
-  virtual Double_t   GetParameter(const char* name) const;
-
-  // read configuration from a file
-  void InitInputs(const char* filename=""){
-    const std::string sf(filename);
-    fInputParameters->configureFromFile(sf);
-  }
-
-  void PrintInputs(ostream& out) const {
-    fInputParameters->print(out);
-  }
-
-  Bool_t InitStarLight() {
-    if (not fInputParameters->init()) {
-      Fatal("InitStarLight", "parameter initialization has failed");
-      fErrorStatus = -1;
-      return false;
-    }
-    return fStarLight->init();
-  }
-
-  void SetParameter(const char* line);
-
-  void SetInput(inputParameters *in) { fInputParameters = in; }
-  inputParameters* GetInput() { return fInputParameters; }
-
-  Int_t GetErrorStatusFlag() const { return fErrorStatus; }
-
-  // boost event to the experiment CM frame
-  void BoostEvent() {
-    fEvent.boost(0.5*(TMath::ACosH(fInputParameters->beam1LorentzGamma()) -
-                     TMath::ACosH(fInputParameters->beam2LorentzGamma())));
-  }
-
- private:
-  bool Stable(Int_t pdgCode) const {
-    switch(TMath::Abs(pdgCode)) {
-    case 11: // electon
-    case 12: // e_neutreno
-    case 13: // Muon
-    case 14: // Muon_neutreno
-    case 16: // tau_neutreno
-    case 22: // Photon
-    case 211:// Charge Pion
-    case 130:// K0 long
-    case 310:// K0 short
-    case 311:// K0
-    case 321:// Charged K
-    case 2212:// Proton
-    case 2112:// neutron
-    case 3122:// Lamda
-    case 3222:// Sigma+
-    case 3112:// Sigma-
-    case 3322:// Exi0
-    case 3312:// Exi-
-    case 3334:// Omega-
-      return kTRUE; // Stable enough.
-    default:
-      return kFALSE; // Not Stable.
-    } // end switch
-    return kFALSE;
-  }
-  Int_t            fErrorStatus;      //  Error status flag. 0=OK
-  TString          fConfigFileName;   //  Input Configuration file name
-  starlight       *fStarLight;        //! Simulation Class
-  inputParameters *fInputParameters;  //! Pointer to simulation input information.
-  upcEvent         fEvent;            //! object holding STARlight simulated event.
-
-  ClassDef(TStarLight,1); // STARlight interface to ROOT's Virtual Monte Carlo
-} ;
-
-#endif
-
diff --git a/STARLIGHT/starlight/TStarLight/testsl.C~ b/STARLIGHT/starlight/TStarLight/testsl.C~
deleted file mode 100644 (file)
index 2a0becd..0000000
+++ /dev/null
@@ -1,58 +0,0 @@
-// -*- C++ -*-
-// $Id$
-
-void testsl() {
-  Printf("Root.UseThreads=%d", gEnv->GetValue("Root.UseThreads", -1));
-
-  gSystem->Load("libStarLight");
-  gSystem->Load("libTStarLight.so");
-
-  TStarLight* sl = new TStarLight("starlight generator", "title", "");
-
-  sl->SetParameter("BEAM_1_Z     =   82    #Z of projectile");
-  sl->SetParameter("BEAM_1_A     =  208    #A of projectile");
-  sl->SetParameter("BEAM_2_Z     =   82    #Z of target");
-  sl->SetParameter("BEAM_2_A     =  208    #A of target");
-  sl->SetParameter("BEAM_1_GAMMA = 1470    #Gamma of the colliding ions");
-  sl->SetParameter("BEAM_2_GAMMA = 1470    #Gamma of the colliding ions");
-  sl->SetParameter("W_MAX        =   12.0  #Max value of w");
-  sl->SetParameter("W_MIN        =    2.0  #Min value of w");
-  sl->SetParameter("W_N_BINS     =   40    #Bins i w");
-  sl->SetParameter("RAP_MAX      =    8.   #max y");
-  sl->SetParameter("RAP_N_BINS   =   80    #Bins i y");
-  sl->SetParameter("CUT_PT       =    0    #Cut in pT? 0 = (no, 1 = yes)");
-  sl->SetParameter("PT_MIN       =    1.0  #Minimum pT in GeV");
-  sl->SetParameter("PT_MAX       =    3.0  #Maximum pT in GeV");
-  sl->SetParameter("CUT_ETA      =    0    #Cut in pseudorapidity? (0 = no, 1 = yes)");
-  sl->SetParameter("ETA_MIN      =  -10    #Minimum pseudorapidity");
-  sl->SetParameter("ETA_MAX      =   10    #Maximum pseudorapidity");
-  sl->SetParameter("PROD_MODE    =    2    #gg or gP switch (1 = 2-photon, 2 = coherent vector meson (narrow), 3 = coherent vector meson (wide), # 4 = incoherent vector meson, 5 = A+A DPMJet single, 6 = A+A DPMJet double, 7 = p+A DPMJet single, 8 = p+A Pythia single )");
-  // is N_EVENTS valid
-  sl->SetParameter("N_EVENTS     =   10    #Number of events");
-  sl->SetParameter("PROD_PID     =  113    #Channel of interest (not relevant for photonuclear processes)");
-  sl->SetParameter("RND_SEED     = 34533   #Random number seed");
-  sl->SetParameter("OUTPUT_FORMAT =   2    #Form of the output");
-  sl->SetParameter("BREAKUP_MODE  =   5    #Controls the nuclear breakup");
-  sl->SetParameter("INTERFERENCE  =   0    #Interference (0 = off, 1 = on)");
-  sl->SetParameter("IF_STRENGTH   =   1.   #% of intefernce (0.0 - 0.1)");
-  sl->SetParameter("COHERENT      =   1    #Coherent=1,Incoherent=0");
-  sl->SetParameter("INCO_FACTOR   =   1.   #percentage of incoherence");
-  sl->SetParameter("BFORD         =   9.5  #");
-  sl->SetParameter("INT_PT_MAX    =   0.24 #Maximum pt considered, when interference is turned on");
-  sl->SetParameter("INT_PT_N_BINS = 120    #Number of pt bins when interference is turned on");
-
-  sl->InitStarLight();
-  sl->PrintInputs(std::cout);
-
-  sl->GetParameter("INT_PT_MAX");
-
-  TClonesArray tca("TParticle", 1000);
-
-  for (Int_t counter(0); counter<10; ++counter) {
-    Printf("--------------------------------------------------------------------------------");
-    sl->GenerateEvent();    
-//     sl->ImportParticles(&tca, "ALL");
-//      for (Int_t i=0; i<tca.GetEntries(); ++i)
-//       tca.At(i)->Print();
-  }
-}
diff --git a/STARLIGHT/starlight/include/inputParameters.h~ b/STARLIGHT/starlight/include/inputParameters.h~
deleted file mode 100644 (file)
index 6b06fc7..0000000
+++ /dev/null
@@ -1,412 +0,0 @@
-///////////////////////////////////////////////////////////////////////////
-//
-//    Copyright 2010
-//
-//    This file is part of starlight.
-//
-//    starlight 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 3 of the License, or
-//    (at your option) any later version.
-//
-//    starlight 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 starlight. If not, see <http://www.gnu.org/licenses/>.
-//
-///////////////////////////////////////////////////////////////////////////
-//
-// File and Version Information:
-// $Rev:: 163                         $: revision of last commit
-// $Author:: odjuvsla                 $: author of last commit
-// $Date:: 2013-10-06 16:18:06 +0200 #$: date of last commit
-//
-// Description:
-//
-//
-//
-///////////////////////////////////////////////////////////////////////////
-
-
-#ifndef INPUTPARAMETERS_H
-#define INPUTPARAMETERS_H
-
-
-#include "starlightconstants.h"
-#include "inputParser.h"
-#include "singleton.h"
-#include <string>
-#include <ostream>
-#include <vector>
-#include <sstream>
-
-class parameterbase;
-
-
-class parameterlist
-{
-public:
-
-    parameterlist() : _parameters(0) {}
-
-    void add(parameterbase* p) {
-        _parameters.push_back(p);
-    }
-
-    // Returns a string with a key of the current state of the parameter list
-    // only
-    inline std::string validationKey();
-    
-
-private:
-
-    std::vector<parameterbase*> _parameters;
-
-};
-
-// Base class for parameters, needed to keep a list of parameters
-class parameterbase
-{
-public:
-
-    // Add this to parameter list
-    parameterbase()
-    {
-        _parameters.add(this);
-    }
-    virtual ~parameterbase() {}
-
-    virtual std::string validationkey() = 0;
-
-    template<typename T>
-    std::string toString(T v)
-    {
-        std::stringstream s;
-        s << v;
-        return s.str();
-    }
-    inline friend std::ostream& operator<<(std::ostream& os, const parameterbase& par);
-    // List of all parameters
-    static parameterlist _parameters;
-
-
-   
-};
-// Need to init the static variable
-// parameterlist parameterbase::_parameters;
-
-
-// The actual parameter class
-// validate parameter specifies if the parameter should be a part of the validity check of the current parameters
-template<typename T, bool validate>
-class parameter : public parameterbase
-{
-public:
-
-    // Constructor
-    parameter(const std::string &name, T value, bool required = true) :parameterbase(),_name(name), _value(value), _validate(validate), _required(required) {}
-
-    virtual ~parameter() {}
-//     T operator()() const {
-//         return _value;
-//     }
-
-    parameter &operator=(T v) { _value = v; return *this;}
-    T* ptr() const {
-        return const_cast<T*>(&_value);
-    }
-    
-    T value() const { return _value; }
-    
-    std::string name() const { return _name;}
-    
-    bool required() const { return _required; }
-    
-    void setValue(T v) { _value = v; }
-    
-    void setName(std::string name) { _name = name; }
-    
-    void setRequired(bool r) { _required = r; }
-    
-    // Validation key for this parameter
-    std::string validationkey()
-    {
-        return (_validate ? _name + ":" + toString(_value) + "-" : std::string(""));
-    }
-
-    template<typename S, bool v>
-    inline friend std::ostream& operator<<(std::ostream& os, const parameter<S,v>& par);
-
-
-
-private:
-    std::string _name;
-
-    T _value; // Value
-    bool _validate; // true if a change in the parameter invalidates x-sec tables
-    bool _required; // true if this is required option.
-
-    parameter();
-};
-
-template<typename S, bool v>
-std::ostream& operator<<(std::ostream& os, const parameter<S,v>& par)
-{
-    os << par._value;
-    return os;
-}
-
-std::ostream& operator<<(std::ostream& os, const parameterbase& par)
-{
-    os << par._parameters.validationKey(); 
-    return os;
-}
-std::string parameterlist::validationKey()
-{
-    std::stringstream s;
-    for(unsigned int i = 0; i < _parameters.size(); ++i)
-    {
-        s << _parameters[i]->validationkey(); // Will print names and values of validation parameters
-    }
-    return s.str();
-}
-
-class inputParameters {
-
-private:
-       // inputParameters is now a singleton
-       friend class Singleton<inputParameters>;
-       inputParameters();
-public:
-
-       ~inputParameters();
-
-       bool init();
-       bool configureFromFile(const std::string &configFileName = "./config/slight.in");
-
-       unsigned int beam1Z                () const { return _beam1Z.value();                 }  ///< returns atomic number of beam particle 1
-       unsigned int beam1A                () const { return _beam1A.value();                 }  ///< returns atomic mass number of beam particle 1
-       unsigned int beam2Z                () const { return _beam2Z.value();                 }  ///< returns atomic number of beam particle 2
-       unsigned int beam2A                () const { return _beam2A.value();                 }  ///< returns atomic mass number of beam particle 2
-       double       beamLorentzGamma      () const { return _beamLorentzGamma;               }  ///< returns Lorentz gamma factor of both beams in beam CMS frame
-       double       beam1LorentzGamma     () const { return _beam1LorentzGamma.value();      }  ///< returns Lorentz gamma factor of beam 1 in collider frame
-       double       beam2LorentzGamma     () const { return _beam2LorentzGamma.value();      }  ///< returns Lorentz gamma factor of beam 2 in collider frame
-       double       maxW                  () const { return _maxW.value();                   }  ///< returns maximum mass W of produced hadronic system [GeV/c^2]
-       double       minW                  () const { return _minW.value();                   }  ///< returns minimum mass W of produced hadronic system [GeV/c^2]
-       unsigned int nmbWBins              () const { return _nmbWBins.value();               }  ///< returns number of W bins in lookup table
-       double       maxRapidity           () const { return _maxRapidity.value();            }  ///< returns maximum absolute value of rapidity
-       unsigned int nmbRapidityBins       () const { return _nmbRapidityBins.value();        }  ///< returns number of rapidity bins in lookup table
-       bool         ptCutEnabled          () const { return _ptCutEnabled.value();           }  ///< returns cut in pt
-       double       ptCutMin              () const { return _ptCutMin.value();               }  ///< returns minimum pt
-       double       ptCutMax              () const { return _ptCutMax.value();               }  ///< returns maximum pt
-       bool         etaCutEnabled         () const { return _etaCutEnabled.value();          }  ///< returns cut in eta
-       double       etaCutMin             () const { return _etaCutMin.value();              }  ///< returns minimum eta
-       double       etaCutMax             () const { return _etaCutMax.value();              }  ///< returns maximum eta
-       int          productionMode        () const { return _productionMode.value();         }  ///< returns production mode
-       unsigned int nmbEvents             () const { return _nmbEventsTot.value();           }  ///< returns total number of events to generate
-       int          prodParticleId        () const { return _prodParticleId.value();         }  ///< returns PDG particle ID of produced particle
-       int          randomSeed            () const { return _randomSeed.value();             }  ///< returns seed for random number generator
-       int          outputFormat          () const { return _outputFormat.value();           }  ///< returns output format
-       int          beamBreakupMode       () const { return _beamBreakupMode.value();        }  ///< returns breakup mode for beam particles
-       bool         interferenceEnabled   () const { return _interferenceEnabled.value();    }  ///< returns whether interference is taken into account
-       double       interferenceStrength  () const { return _interferenceStrength.value();   }  ///< returns percentage of interference
-       bool         coherentProduction    () const { return _coherentProduction.value();     }  ///< returns whether production is coherent or incoherent
-       double       incoherentFactor      () const { return _incoherentFactor.value();       }  ///< returns incoherent contribution in vector meson production
-       double       deuteronSlopePar      () const { return _deuteronSlopePar.value();       }  ///< returns slope parameter for deuteron form factor [(GeV/c)^{-2}]
-       double       maxPtInterference     () const { return _maxPtInterference.value();      }  ///< returns maximum p_T for interference calculation [GeV/c]
-       int          nmbPtBinsInterference () const { return _nmbPtBinsInterference.value();  }  ///< returns number of p_T bins for interference calculation
-       double       ptBinWidthInterference() const { return _ptBinWidthInterference.value(); }  ///< returns width of p_T bins for interference calculation [GeV/c]
-       double       minGammaEnergy        () const { return _minGammaEnergy.value();         }  ///< returns minimum gamma energy in case of photo nuclear processes [GeV]
-       double       maxGammaEnergy        () const { return _maxGammaEnergy.value();         }  ///< returns maximum gamma energy in case of photo nuclear processes [GeV]
-       std::string  pythiaParams          () const { return _pythiaParams.value();           }  ///< returns parameters to be passed to pythia
-       bool         pythiaFullEventRecord () const { return _pythiaFullEventRecord.value();  }  ///< returns if the full pythia event record should be printed
-       int          xsecCalcMethod        () const { return _xsecCalcMethod.value();         }  ///< returns the method used for the x-sec calculation
-       int          nThreads              () const { return _nThreads.value();               }  ///< returns the number of threads in case method 1 is used for the x-sec calc
-       unsigned int nBinsQKniehl          () const { return _nBinsQKniehl.value();           }  ///< Number of bins in Q used for the transformation to the impact paramter space of the Kniehl function
-        unsigned int nBinsEKniehl          () const { return _nBinsEKniehl.value();           }  ///< Number of bins in photon energy used for the Kniehl function
-        unsigned int nBinsBKniehl          () const { return _nBinsBKniehl.value();           }  ///< Number of bins in impact parameter used for the Kniehl function
-        double       qMaxKniehl            () const { return _qMaxKniehl.value();             }  ///< Max value of Q used for the Kniehl funcion
-        double       eGammaMinKniehl       () const { return _eGammaMinKniehl.value();        }  ///< Min value of gamma energy used for the Kniehl funcion
-        double       eGammaMaxKniehl       () const { return _eGammaMaxKniehl.value();        }  ///< Max value of gamma energy used for the Kniehl funcion
-        double       bMinKniehl            () const { return _bMinKniehl.value();             }  ///< Min value of impact parameter used for the Kniehl funcion
-        double       bMaxKniehl            () const { return _bMaxKniehl.value();             }  ///< Max value of impact parameter used for the Kniehl funcion
-        
-       starlightConstants::particleTypeEnum    prodParticleType     () const { return _particleType;    }  ///< returns type of produced particle
-       starlightConstants::decayTypeEnum       prodParticleDecayType() const { return _decayType;       }  ///< returns decay type of produced particle
-       starlightConstants::interactionTypeEnum interactionType      () const { return _interactionType; }  ///< returns interaction type
-       // double vmPhotonCoupling();
-       // double slopeParameter();
-       double protonEnergy                () const { return _protonEnergy.value(); }
-
-       void setBeam1Z                (unsigned int v)  {  _beam1Z = v;                 }  ///< returns atomic number of beam particle 1
-       void setBeam1A                (unsigned int v)  {  _beam1A = v;                 }  ///< returns atomic mass number of beam particle 1
-       void setBeam2Z                (unsigned int v)  {  _beam2Z = v;                 }  ///< returns atomic number of beam particle 2
-       void setBeam2A                (unsigned int v)  {  _beam2A = v;                 }  ///< returns atomic mass number of beam particle 2
-       void setBeamLorentzGamma      (double v)  {  _beamLorentzGamma = v;       }  ///< returns Lorentz gamma factor of both beams in beam CMS frame
-       void setBeam1LorentzGamma     (double v)  {  _beam1LorentzGamma = v;      }  ///< returns Lorentz gamma factor of beam 1 in collider frame
-       void setBeam2LorentzGamma     (double v)  {  _beam2LorentzGamma = v;      }  ///< returns Lorentz gamma factor of beam 2 in collider frame
-       void setMaxW                  (double v)  {  _maxW = v;                   }  ///< returns maximum mass W of produced hadronic system [GeV/c^2]
-       void setMinW                  (double v)  {  _minW = v;                   }  ///< returns minimum mass W of produced hadronic system [GeV/c^2]
-       void setNmbWBins              (unsigned int v)  {  _nmbWBins = v;               }  ///< returns number of W bins in lookup table
-       void setMaxRapidity           (double v)  {  _maxRapidity = v;            }  ///< returns maximum absolute value of rapidity
-       void setNmbRapidityBins       (unsigned int v)  {  _nmbRapidityBins = v;        }  ///< returns number of rapidity bins in lookup table
-       void setPtCutEnabled          (bool v)  {  _ptCutEnabled = v;           }  ///< returns cut in pt
-       void setPtCutMin              (double v)  {  _ptCutMin = v;               }  ///< returns minimum pt
-       void setPtCutMax              (double v)  {  _ptCutMax = v;               }  ///< returns maximum pt
-       void setEtaCutEnabled         (bool v)  {  _etaCutEnabled = v;          }  ///< returns cut in eta
-       void setEtaCutMin             (double v)  {  _etaCutMin = v;              }  ///< returns minimum eta
-       void setEtaCutMax             (double v)  {  _etaCutMax = v;              }  ///< returns maximum eta
-       void setProductionMode        (int v)  {  _productionMode = v;         }  ///< returns production mode
-       void setNmbEvents             (unsigned int v)  {  _nmbEventsTot = v;           }  ///< returns total number of events to generate
-       void setProdParticleId        (int v)  {  _prodParticleId = v;         }  ///< returns PDG particle ID of produced particle
-       void setRandomSeed            (int v)  {  _randomSeed = v;             }  ///< returns seed for random number generator
-       void setOutputFormat          (int v)  {  _outputFormat = v;           }  ///< returns output format
-       void setBeamBreakupMode       (int v)  {  _beamBreakupMode = v;        }  ///< returns breakup mode for beam particles
-       void setInterferenceEnabled   (bool v)  {  _interferenceEnabled = v;    }  ///< returns whether interference is taken into account
-       void setInterferenceStrength  (double v)  {  _interferenceStrength = v;   }  ///< returns percentage of interference
-       void setCoherentProduction    (bool v)  {  _coherentProduction = v;     }  ///< returns whether production is coherent or incoherent
-       void setIncoherentFactor      (double v)  {  _incoherentFactor = v;       }  ///< returns incoherent contribution in vector meson production
-       void setDeuteronSlopePar      (double v)  {  _deuteronSlopePar = v;       }  ///< returns slope parameter for deuteron form factor [(GeV/c)^{-2}]
-       void setMaxPtInterference     (double v)  {  _maxPtInterference = v;      }  ///< returns maximum p_T for voiderference calculation [GeV/c]
-       void setNmbPtBinsInterference (int v)  {  _nmbPtBinsInterference = v;  }  ///< returns number of p_T bins for interference calculation
-       void setPtBinWidthInterference(double v)  {  _ptBinWidthInterference = v; }  ///< returns width of p_T bins for voiderference calculation [GeV/c]
-       void setMinGammaEnergy        (double v)  {  _minGammaEnergy = v;         }  ///< returns minimum gamma energy in case of photo nuclear processes [GeV]
-       void setMaxGammaEnergy        (double v)  {  _maxGammaEnergy = v;         }  ///< returns maximum gamma energy in case of photo nuclear processes [GeV]
-       void setPythiaParams          (std::string v)  {  _pythiaParams = v;           }  ///< returns parameters to be passed to pythia
-       void setPythiaFullEventRecord (bool v)  {  _pythiaFullEventRecord = v;  }  ///< returns if the full pythia event record should be prvoided
-       void setXsecCalcMethod        (int v)  {  _xsecCalcMethod = v;         }  ///< returns the method used for the x-sec calculation
-       void setNThreads              (int v)  {  _nThreads = v;               }  ///< returns the number of threads in case method 1 is used for the x-sec calc
-       void setNBinsQKniehl          (unsigned int v)  {  _nBinsQKniehl = v;           }  ///< Number of bins in Q used for the transformation to the impact paramter space of the Kniehl function
-        void setNBinsEKniehl          (unsigned int v)  {  _nBinsEKniehl = v;           }  ///< Number of bins in photon energy used for the Kniehl function
-        void setNBinsBKniehl          (unsigned int v)  {  _nBinsBKniehl = v;           }  ///< Number of bins in impact parameter used for the Kniehl function
-        void setQMaxKniehl            (double v)  {  _qMaxKniehl = v;             }  ///< Max value of Q used for the Kniehl funcion
-        void setEGammaMinKniehl       (double v)  {  _eGammaMinKniehl = v;        }  ///< Min value of gamma energy used for the Kniehl funcion
-        void setEGammaMaxKniehl       (double v)  {  _eGammaMaxKniehl = v;        }  ///< Max value of gamma energy used for the Kniehl funcion
-        void setBMinKniehl            (double v)  {  _bMinKniehl = v;             }  ///< Min value of impact parameter used for the Kniehl funcion
-        void setBMaxKniehl            (double v)  {  _bMaxKniehl = v;             }  ///< Max value of impact parameter used for the Kniehl funcion
-        
-       void setProdParticleType      (starlightConstants::particleTypeEnum v)   { _particleType = v;    }  ///< returns type of produced particle
-       void setProdParticleDecayType (starlightConstants::decayTypeEnum v)        { _decayType = v;       }  ///< returns decay type of produced particle
-       void setInteractionType       (starlightConstants::interactionTypeEnum v)  { _interactionType = v; }  ///< returns interaction type
-        
-       // double vmPhotonCoupling();
-       // double slopeParameter();
-       void setProtonEnergy        (double v)  { _protonEnergy = v; }
-       
-//     template<typename T>
-       inline bool setParameter(std::string expression);
-       
-       std::ostream& print(std::ostream& out) const;  ///< prints parameter summary
-       std::ostream& write(std::ostream& out) const;  ///< writes parameters back to an ostream
-       
-       std::string parameterValueKey() const; ///< Generates key for the current parameters
-
-  
-private:
-
-    
-// To indicate if the crossection table should be re-calculated if parameter changes
-#define VALIDITY_CHECK true
-#define NO_VALIDITY_CHECK false
-       
-       std::string _configFileName;  ///< path to configuration file (default = ./config/slight.in)
-
-       // config file parameters
-       parameter<unsigned int,VALIDITY_CHECK>     _beam1Z;                  ///< atomic number of beam particle 1
-       parameter<unsigned int,VALIDITY_CHECK>     _beam1A;                  ///< atomic mass number of beam particle 1
-       parameter<unsigned int,VALIDITY_CHECK>     _beam2Z;                  ///< atomic number of beam particle 2
-       parameter<unsigned int,VALIDITY_CHECK>     _beam2A;                  ///< atomic mass number of beam particle 2
-       parameter<double, VALIDITY_CHECK>          _beam1LorentzGamma;       ///< Lorentz gamma factor of beam 1 in collider frame
-       parameter<double, VALIDITY_CHECK>          _beam2LorentzGamma;       ///< Lorentz gamma factor of beam 2 in collider frame
-       parameter<double, VALIDITY_CHECK>          _maxW;                    ///< maximum mass W of produced hadronic system [GeV/c^2]
-       parameter<double, VALIDITY_CHECK>          _minW;                    ///< minimum mass W of produced hadronic system; if set to -1 default value is taken [GeV/c^2]
-       parameter<unsigned int, VALIDITY_CHECK>    _nmbWBins;                ///< number of W bins in lookup table
-       parameter<double, VALIDITY_CHECK>          _maxRapidity;             ///< maximum absolute value of rapidity
-       parameter<unsigned int, VALIDITY_CHECK>    _nmbRapidityBins;         ///< number of rapidity bins in lookup table
-       parameter<bool, VALIDITY_CHECK>            _ptCutEnabled;            ///< en/disables cut in pt
-       parameter<double, VALIDITY_CHECK>          _ptCutMin;                ///< minimum pt, if cut is enabled
-       parameter<double, VALIDITY_CHECK>          _ptCutMax;                ///< maximum pt, if cut is enabled
-       parameter<bool, VALIDITY_CHECK>            _etaCutEnabled;           ///< en/disables cut in eta
-       parameter<double, VALIDITY_CHECK>          _etaCutMin;               ///< minimum eta, if cut is enabled
-       parameter<double, VALIDITY_CHECK>          _etaCutMax;               ///< maximum eta, if cut is enabled
-       parameter<unsigned int, VALIDITY_CHECK>    _productionMode;          ///< \brief production mode
-                                                                            ///<
-                                                                            ///< 1 = photon-photon fusion,
-                                                                            ///< 2 = narrow vector meson resonance in photon-Pomeron fusion,
-                                                                            ///< 3 = Breit-Wigner vector meson resonance in photon-Pomeron fusion
-       parameter<unsigned int, VALIDITY_CHECK>    _nmbEventsTot;            ///< total number of events to generate
-       parameter<unsigned int, VALIDITY_CHECK>    _prodParticleId;          ///< PDG particle ID of produced particle
-       parameter<unsigned int, VALIDITY_CHECK>    _randomSeed;              ///< seed for random number generator
-       parameter<unsigned int, NO_VALIDITY_CHECK> _outputFormat;            ///< \brief output format
-                                                                            ///<
-                                                                            ///< 1 = ASCII
-                                                                            ///< 2 = GSTARtext,
-                                                                            ///< 3 = PAW ntuple (not working)
-       parameter<unsigned int, VALIDITY_CHECK>    _beamBreakupMode;         ///< \brief breakup mode for beam particles
-                                                                            ///<
-                                                                            ///< 1 = hard sphere nuclei (b > 2R),
-                                                                            ///< 2 = both nuclei break up (XnXn),
-                                                                            ///< 3 = a single neutron from each nucleus (1n1n),
-                                                                            ///< 4 = neither nucleon breaks up (with b > 2R),
-                                                                            ///< 5 = no hadronic break up (similar to option 1, but with the actual hadronic interaction)
-       parameter<bool, VALIDITY_CHECK>            _interferenceEnabled;     ///< if VALIDITY_CHECK, interference is taken into account
-       parameter<double, VALIDITY_CHECK>          _interferenceStrength;    ///< percentage of interference: from 0 = none to 1 = full
-       parameter<bool, VALIDITY_CHECK>            _coherentProduction;      ///< if VALIDITY_CHECK, production is coherent, else incoherent
-       parameter<double, VALIDITY_CHECK>          _incoherentFactor;        ///< allows to scale the incoherent contribution in vector meson production
-       parameter<double, VALIDITY_CHECK>          _deuteronSlopePar;        ///< slope parameter for deuteron form factor [(GeV/c)^{-2}]
-       parameter<double, VALIDITY_CHECK>          _maxPtInterference;       ///< maximum p_T for interference calculation [GeV/c]
-       parameter<unsigned int, VALIDITY_CHECK>    _nmbPtBinsInterference;   ///< number of p_T bins for interference calculation
-       parameter<double, VALIDITY_CHECK>          _ptBinWidthInterference;  ///< width of p_T bins for interference calculation [GeV/c]
-       parameter<double, VALIDITY_CHECK>          _protonEnergy;
-       parameter<double, VALIDITY_CHECK>          _minGammaEnergy;          ///< minimum gamma energy in case of photo nuclear processes [GeV]
-       parameter<double, VALIDITY_CHECK>          _maxGammaEnergy;          ///< maximum gamma energy in case of photo nuclear processes [GeV]
-       parameter<std::string,NO_VALIDITY_CHECK>   _pythiaParams;            ///< semi-colon separated parameters to pass to pythia, e.g. "mstj(1)=0;paru(13)=0.1" 
-       parameter<bool, NO_VALIDITY_CHECK>         _pythiaFullEventRecord;   ///< if the full pythia event record should be in the outputu
-       parameter<unsigned int, VALIDITY_CHECK>    _xsecCalcMethod;          ///< Select x-sec calc method. (0 is standard starlight method, 1 must be used for assym. collisions (e.g. p-A), but is slow)
-       parameter<unsigned int, NO_VALIDITY_CHECK> _nThreads;                ///< Number of threads used in the case of using method 1 for calculating the x-sections
-        parameter<unsigned int, VALIDITY_CHECK>    _nBinsQKniehl;            ///< Number of bins in Q used for the transformation to the impact paramter space of the Kniehl function
-        parameter<unsigned int, VALIDITY_CHECK>    _nBinsEKniehl;            ///< Number of bins in photon energy used for the Kniehl function
-        parameter<unsigned int, VALIDITY_CHECK>    _nBinsBKniehl;            ///< Number of bins in impact parameter used for the Kniehl function
-        parameter<double, VALIDITY_CHECK>          _qMaxKniehl;              ///< Max value of Q used for the Kniehl funcion
-        parameter<double, VALIDITY_CHECK>          _eGammaMinKniehl;         ///< Min value of gamma energy used for the Kniehl funcion
-        parameter<double, VALIDITY_CHECK>          _eGammaMaxKniehl;         ///< Max value of gamma energy used for the Kniehl funcion
-        parameter<double, VALIDITY_CHECK>          _bMinKniehl;              ///< Min value of impact parameter used for the Kniehl funcion
-        parameter<double, VALIDITY_CHECK>          _bMaxKniehl;              ///< Max value of impact parameter used for the Kniehl funcion
-        
-       
-       starlightConstants::particleTypeEnum       _particleType;
-       starlightConstants::decayTypeEnum          _decayType;
-       starlightConstants::interactionTypeEnum    _interactionType;
-
-       double                         _beamLorentzGamma;                    ///< Lorentz gamma factor of the beams in CMS frame, not an input parameter
-       
-       inputParser _ip;
-       
-};
-
-#define inputParametersInstance Singleton<inputParameters>::instance()
-
-// template<typename T>
-inline 
-bool inputParameters::setParameter(std::string expression)
-{
-   
-    return _ip.parseString(expression);
-   
-   
-}
-
-inline
-std::ostream&
-operator <<(std::ostream&          out,
-            const inputParameters& par)
-{
-       return par.print(out);
-}
-
-#endif  // INPUTPARAMETERS_H
diff --git a/STARLIGHT/starlight/include/inputParser.h~ b/STARLIGHT/starlight/include/inputParser.h~
deleted file mode 100644 (file)
index 174e164..0000000
+++ /dev/null
@@ -1,148 +0,0 @@
-///////////////////////////////////////////////////////////////////////////
-//
-//    Copyright 2011
-//
-//    This file is part of starlight.
-//
-//    starlight 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 3 of the License, or
-//    (at your option) any later version.
-//
-//    starlight 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 starlight. If not, see <http://www.gnu.org/licenses/>.
-//
-///////////////////////////////////////////////////////////////////////////
-//
-// File and Version Information:
-// $Rev:: 155                       $: revision of last commit
-// $Author:: odjuvsla           $: author of last commit
-// $Date:: 2013-10-06 16:17:50 +#$: date of last commit
-//
-// Description:
-//
-//
-//
-///////////////////////////////////////////////////////////////////////////
-
-
-#ifndef INPUTPARSER_H
-#define INPUTPARSER_H
-
-#include <string>
-#include <typeinfo>
-#include <iostream>
-#include <map>
-
-#include <reportingUtils.h>
-
-class inputParser
-{
-public:
-  
-  /** Constructor */
-  inputParser();
-
-  /** Destructor */
-  ~inputParser();
-
-  /** Parse a file */
-  int parseFile(std::string filename);
-
-  /** Parse a file */
-  int parseString(std::string str);
-  
-  /** Add parameter to pass */
-  void addIntParameter(std::string name, int *var, bool required = true);
-
-  /** Add parameter to pass */
-  void addUintParameter(std::string name, unsigned int *var, bool required = true);
-
-  /** Add parameter to pass */
-  void addFloatParameter(std::string name, float *var, bool required = true);
-
-  /** Add parameter to pass */
-  void addDoubleParameter(std::string name, double *var, bool required = true);
-
-  /** Add parameter to pass */
-  void addBoolParameter(std::string name, bool *var, bool required = true);
-  /** Add parameter to pass */
-  void addStringParameter(std::string name, std::string *var, bool required = true);
-  
-  /** Print info */
-  void printParameterInfo(std::ostream &out = std::cout);
-  
-  /** Validate */
-  bool validateParameters(std::ostream &warnOut = std::cout, std::ostream &errOut = std::cerr);
-  
-  /** Add a parameter */
-  template<typename S>
-  inline void addParameter(S &param);
-  
-  /** Add a parameter */
-  template<typename P>
-  inline void addParameter(const std::string &name, P *varPtr, bool required = false);
-
-private:
-  
-  template <class T>
-  class _parameter
-  {
-  public:
-    _parameter(std::string name, T *val, bool required = true, bool found = false) : _name(name), _val(val), _required(required), _found(found){}
-    
-    bool operator==(const _parameter &rhs) const { return _name == rhs._name; }
-    
-    bool operator<(const _parameter &rhs) const { return _name.c_str()[0] < rhs._name.c_str()[0]; }
-    
-    void printParameterInfo(std::ostream &out = std::cout) 
-    {
-      out << std::boolalpha << _name << "\t\t";
-      if(_found)
-      {
-       out << *_val << std::endl;
-      }
-      else
-      {
-       out << "NOT FOUND" << std::endl;
-      }
-      out << std::noboolalpha;
-    }
-    
-    
-    std::string _name;
-    T *_val;
-    bool _required;
-    bool _found;
-  };
-  
-  std::map<std::string, _parameter<int> > _intParameters;
-  std::map<std::string, _parameter<unsigned int> > _uintParameters;
-  std::map<std::string, _parameter<float> > _floatParameters;
-  std::map<std::string, _parameter<double> > _doubleParameters;
-  std::map<std::string, _parameter<bool> > _boolParameters;
-  std::map<std::string, _parameter<std::string> > _stringParameters;
-  
-};
-
-template<typename S>
-void inputParser::addParameter(S& param)
-{
-  addParameter(param.name(), param.ptr(), param.required());
-
-}
-
-template<typename P>
-void inputParser::addParameter(const std::string& name, P* /*varPtr*/, bool /*required*/)
-{
-  printWarn << "Trying to add unknown parameter type with name: " << name;
-}
-
-
-#endif // INPUTPARSER_H
diff --git a/STARLIGHT/starlight/include/nBodyPhaseSpaceGen.h~ b/STARLIGHT/starlight/include/nBodyPhaseSpaceGen.h~
deleted file mode 100644 (file)
index a39446e..0000000
+++ /dev/null
@@ -1,229 +0,0 @@
-///////////////////////////////////////////////////////////////////////////
-//
-//    Copyright 2010
-//
-//    This file is part of starlight.
-//
-//    starlight 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 3 of the License, or
-//    (at your option) any later version.
-//       
-//    starlight 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 starlight. If not, see <http://www.gnu.org/licenses/>.
-//
-///////////////////////////////////////////////////////////////////////////
-//
-// File and Version Information:
-// $Rev:: 157                         $: revision of last commit
-// $Author:: odjuvsla                 $: author of last commit
-// $Date:: 2013-10-06 16:17:54 +0200 #$: date of last commit
-//
-// Description:
-//    calculates n-body phase space (constant matrix element) using various algorithms
-//    
-//    the n-body decay is split up into (n - 2) successive 2-body decays
-//    each 2-body decay is considered in its own center-of-mass frame thereby
-//    separating the mass from the (trivial) angular dependence
-//    
-//    the event is boosted into the same frame in which the n-body system is
-//    given
-//    
-//    based on:
-//    GENBOD (CERNLIB W515), see F. James, "Monte Carlo Phase Space", CERN 68-15 (1968)
-//    NUPHAZ, see M. M. Block, "Monte Carlo phase space evaluation", Comp. Phys. Commun. 69, 459 (1992)
-//    S. U. Chung, "Spin Formalism", CERN Yellow Report
-//    S. U. Chung et. al., "Diffractive Dissociation for COMPASS"
-//    
-//    index convention:
-//    - all vectors have the same size (= number of decay daughters)
-//    - index i corresponds to the respective value in the (i + 1)-body system: effective mass M, break-up momentum, angles
-//    - thus some vector elements are not used like breakupMom[0], theta[0], phi[0], ...
-//      this overhead is negligible compared to the ease of notation
-//    
-//    the following graph illustrates how the n-body decay is decomposed into a sequence of two-body decays
-//    
-//    n-body       ...                   3-body                 2-body                 single daughter
-//    
-//    m[n - 1]                           m[2]                   m[1]
-//     ^                                  ^                      ^
-//     |                                  |                      |
-//     |                                  |                      |
-//    M[n - 1] --> ... -->               M[2] -->               M[1] -->               M    [0] = m[0]
-//    theta[n - 1] ...                   theta[2]               theta[1]               theta[0] = 0 (not used)
-//    phi  [n - 1] ...                   phi  [2]               phi  [1]               phi  [0] = 0 (not used)
-//    mSum [n - 1] ...                   mSum [2]               mSum [1]               mSum [0] = m[0]
-//    = sum_0^(n - 1) m[i]               = m[2] + m[1] + m[0]   = m[1] + m[0]
-//    breakUpMom[n - 1] ...              breakUpMom[2]          breakUpMom[1]          breakUpMom[0] = 0 (not used)
-//    = q(M[n - 1], m[n - 1], M[n - 2])  = q(M[2], m[2], M[1])  = q(M[1], m[1], m[0])
-//
-//
-///////////////////////////////////////////////////////////////////////////
-
-
-#ifndef NBODYPHASESPACEGEN_H
-#define NBODYPHASESPACEGEN_H
-
-
-#include <iostream>
-#include <vector>
-
-#include "reportingUtils.h"
-#include "lorentzvector.h"
-#include "randomgenerator.h"
-#include "starlightconstants.h"
-
-
-// small helper functions
-// calculates factorial
-inline
-unsigned int
-factorial(const unsigned int n)
-{
-       unsigned int fac = 1;
-       for (unsigned int i = 1; i <= n; ++i)
-               fac *= i;
-       return fac;
-}
-
-
-// computes breakup momentum of 2-body decay
-inline
-double
-breakupMomentum(const double M,   // mass of mother particle
-                const double m1,  // mass of daughter particle 1
-                const double m2)  // mass of daughter particle 2
-{
-  if (M < m1 + m2)
-    return 0;
-  return sqrt((M - m1 - m2) * (M + m1 + m2) * (M - m1 + m2) * (M + m1 - m2)) / (2 * M);
-}
-
-
-class nBodyPhaseSpaceGen {
-
-public:
-
-       nBodyPhaseSpaceGen();
-       virtual ~nBodyPhaseSpaceGen();
-  
-       // generator setup
-       /// sets decay constants and prepares internal variables
-       bool setDecay(const std::vector<double>& daughterMasses);  // daughter particle masses
-       bool setDecay(const unsigned int         nmbOfDaughters,   // number of daughter particles
-                     const double*              daughterMasses);  // array of daughter particle masses
-  
-       // random generator
-       double random ()                        { return randyInstance.Rndom(); }  ///< returns number from internal random generator
-
-       // high-level generator interface
-       /// generates full event with certain n-body mass and momentum and returns event weight
-       double generateDecay        (const lorentzVector& nBody);          // Lorentz vector of n-body system in lab frame
-       /// \brief generates full event with certain n-body mass and momentum only when event is accepted (return value = true)
-       /// this function is more efficient, if only weighted events are needed
-       bool   generateDecayAccepted(const lorentzVector& nBody,           // Lorentz vector of n-body system in lab frame
-                                    const double         maxWeight = 0);  // if positive, given value is used as maximum weight, otherwise _maxWeight
-
-       void   setMaxWeight          (const double maxWeight) { _maxWeight = maxWeight;    }  ///< sets maximum weight used for hit-miss MC
-       double maxWeight             () const                 { return _maxWeight;         }  ///< returns maximum weight used for hit-miss MC
-       double normalization         () const                 { return _norm;              }  ///< returns normalization used in weight calculation
-       double eventWeight           () const                 { return _weight;            }  ///< returns weight of generated event
-       double maxWeightObserved     () const                 { return _maxWeightObserved; }  ///< returns maximum observed weight since instantiation
-       void   resetMaxWeightObserved()                       { _maxWeightObserved = 0;    }  ///< sets maximum observed weight back to zero
-
-       /// estimates maximum weight for given n-body mass
-       double estimateMaxWeight(const double       nBodyMass,                 // sic!
-                                const unsigned int nmbOfIterations = 10000);  // number of generated events
-
-       /// \brief applies event weight in form of hit-miss MC
-       /// assumes that event weight has been already calculated by calcWeight()
-       /// if maxWeight > 0 value is used as maximum weight, otherwise _maxWeight value is used
-       inline bool eventAccepted(const double maxWeight = 0);
-
-       //----------------------------------------------------------------------------
-       // trivial accessors
-       const lorentzVector&              daughter        (const int index) const { return _daughters[index];  }  ///< returns Lorentz vector of daughter at index
-       const std::vector<lorentzVector>& daughters       ()                const { return _daughters;         }  ///< returns Lorentz vectors of all daughters
-       unsigned int                      nmbOfDaughters  ()                const { return _n;                 }  ///< returns number of daughters
-       double                            daughterMass    (const int index) const { return _m[index];          }  ///< returns invariant mass of daughter at index
-       double                            intermediateMass(const int index) const { return _M[index];          }  ///< returns intermediate mass of (index + 1)-body system
-       double                            breakupMom      (const int index) const { return _breakupMom[index]; }  ///< returns breakup momentum in (index + 1)-body RF
-       double                            cosTheta        (const int index) const { return _cosTheta[index];   }  ///< returns polar angle in (index + 1)-body RF
-       double                            phi             (const int index) const { return _phi[index];        }  ///< returns azimuth in (index + 1)-body RF
-
-
-       std::ostream& print(std::ostream& out = std::cout) const;  ///< prints generator status
-       friend std::ostream& operator << (std::ostream&             out,
-                                         const nBodyPhaseSpaceGen& gen)
-       { return gen.print(out); }
-
-private:
-
-       //----------------------------------------------------------------------------
-       // low-level generator interface
-       /// randomly choses the (n - 2) effective masses of the respective (i + 1)-body systems
-       void pickMasses(const double nBodyMass);  // total energy of n-body system in its RF
-
-       /// \brief computes event weight and breakup momenta
-       /// operates on vector of intermediate two-body masses prepared by pickMasses()
-       double calcWeight();
-
-       /// randomly choses the (n - 1) polar and (n - 1) azimuthal angles in the respective (i + 1)-body RFs
-       inline void pickAngles();
-               
-       /// \brief calculates full event kinematics from the effective masses of the (i + 1)-body systems and the Lorentz vector of the decaying system
-       /// uses the break-up momenta calculated by calcWeight() and angles from pickAngles()
-       void calcEventKinematics(const lorentzVector& nBody);  // Lorentz vector of n-body system in lab frame
-
-       // external parameters
-       std::vector<double> _m;  ///< masses of daughter particles
-
-       // internal variables
-       unsigned int               _n;                  ///< number of daughter particles
-       std::vector<double>        _M;                  ///< effective masses of (i + 1)-body systems
-       std::vector<double>        _cosTheta;           ///< cosine of polar angle of the 2-body decay of the (i + 1)-body system
-       std::vector<double>        _phi;                ///< azimuthal angle of the 2-body decay of the (i + 1)-body system
-       std::vector<double>        _mSum;               ///< sums of daughter particle masses
-       std::vector<double>        _breakupMom;         ///< breakup momenta for the two-body decays: (i + 1)-body --> daughter_(i + 1) + i-body
-       std::vector<lorentzVector> _daughters;          ///< Lorentz vectors of the daughter particles
-       double                     _norm;               ///< normalization value
-       double                     _weight;             ///< phase space weight of generated event
-       double                     _maxWeightObserved;  ///< maximum event weight calculated processing the input data
-       double                     _maxWeight;          ///< maximum weight used to weight events in hit-miss MC
-
-
-};
-
-
-inline
-void
-nBodyPhaseSpaceGen::pickAngles()
-{
-       for (unsigned int i = 1; i < _n; ++i) {  // loop over 2- to n-bodies
-               _cosTheta[i] = 2 * random() - 1;  // range [-1,    1]
-               _phi[i]      = starlightConstants::twoPi * random();  // range [ 0, 2 pi]
-       }
-}
-
-
-inline 
-bool
-nBodyPhaseSpaceGen::eventAccepted(const double maxWeight)  // if maxWeight > 0, given value is used as maximum weight, otherwise _maxWeight
-{
-       const double max = (maxWeight <= 0) ? _maxWeight : maxWeight;
-       if (max <= 0) {
-               printWarn << "maximum weight = " << max << " does not make sense. rejecting event." << std::endl;
-               return false;
-       }
-       if ((_weight / max) > random())
-               return true;
-       return false;
-}
-
-
-#endif  // NBODYPHASESPACEGEN_H
diff --git a/STARLIGHT/starlight/include/singleton.h~ b/STARLIGHT/starlight/include/singleton.h~
deleted file mode 100644 (file)
index 4466374..0000000
+++ /dev/null
@@ -1,73 +0,0 @@
-#ifndef SINGLETON_H
-#define SINGLETON_H
-
-#ifdef CPP11
-#include <atomic>
-#include <mutex>
-
-
-template<class T>
-class Singleton {
-public:
-  static T& instance()
-  {
-    T * tmp = instance_.load(std::memory_order_consume);
-    if (!tmp) {
-      std::lock_guard<std::mutex> guard(instantiation_mutex);
-      tmp = instance_.load(std::memory_order_consume);
-      if (!tmp) {
-        tmp = new T;
-        instance_.store(tmp, std::memory_order_release);
-      }
-    }
-    return *tmp;
-  }
-private:
-  static std::atomic<T *> instance_;
-  static std::mutex instantiation_mutex;
-};
-
-template<class T>
-std::atomic<T *> Singleton<T>::instance_(0);
-
-template<class T>
-std::mutex Singleton<T>::instantiation_mutex;
-
-#else
-
-#include "slmutex.h"
-
-template<class T>
-class Singleton {
-public:
-  static T& instance()
-  {
-    T *tmp = _instance;
-    if (!tmp) {
-      Lockguard<MutexPosix> guard(&_instantiation_mutex);
-      if (!tmp) {
-        tmp = new T;
-       _instance = tmp;
-      }
-    }
-    return *tmp;
-  }
-private:
-  
-  static T * _instance;
-  
-  static MutexPosix _instantiation_mutex;
-  
-};
-
-template<class T>
-T *Singleton<T>::_instance(0);
-
-template<class T>
-MutexPosix Singleton<T>::_instantiation_mutex;
-
-#endif
-
-#endif
-
-  
\ No newline at end of file
diff --git a/STARLIGHT/starlight/include/slmutex.h~ b/STARLIGHT/starlight/include/slmutex.h~
deleted file mode 100644 (file)
index 0035392..0000000
+++ /dev/null
@@ -1,37 +0,0 @@
-#ifndef SLMUTEX_H
-#define SLMUTEX_H
-
-#include <pthread.h>
-
-class MutexPosix
-{
-private:
-  MutexPosix(const MutexPosix &);
-  
-  pthread_mutex_t _mutex;
-  
-public:
-  MutexPosix() {}
-  int lock() { return pthread_mutex_lock(&_mutex); }
-  int unlock() { return pthread_mutex_unlock(&_mutex); }
-  
-};
-
-template<typename M>
-class Lockguard
-{
-private:
-  
-  M *_mutex;
-  Lockguard(const Lockguard & guard); // Do not implement
-public:
-  
-  Lockguard(M *mutex): _mutex(mutex) {_mutex->lock();}
-  ~Lockguard() { _mutex->unlock(); }
-  
-};
-
-#endif
\ No newline at end of file
diff --git a/STARLIGHT/starlight/include/starlightStandalone.h~ b/STARLIGHT/starlight/include/starlightStandalone.h~
deleted file mode 100644 (file)
index 505ad80..0000000
+++ /dev/null
@@ -1,77 +0,0 @@
-///////////////////////////////////////////////////////////////////////////
-//
-//    Copyright 2010
-//
-//    This file is part of starlight.
-//
-//    starlight 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 3 of the License, or
-//    (at your option) any later version.
-//
-//    starlight 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 starlight. If not, see <http://www.gnu.org/licenses/>.
-//
-///////////////////////////////////////////////////////////////////////////
-//
-// File and Version Information:
-// $Rev:: 102                         $: revision of last commit
-// $Author:: odjuvsla                 $: author of last commit
-// $Date:: 2012-10-22 23:25:54 +0200 #$: date of last commit
-//
-// Description:
-//
-//
-//
-///////////////////////////////////////////////////////////////////////////
-
-
-#ifndef STARLIGHTSTANDALONE_H
-#define STARLIGHTSTANDALONE_H
-
-
-#include <string>
-
-class upcEvent;
-
-class starlight;
-class inputParameters;
-
-
-class starlightStandalone {
-
-public:
-
-       starlightStandalone();
-       ~starlightStandalone();
-
-       bool init();  ///< reads configuration file and initializes startlight
-       bool run ();  ///< creates output file and runs starlight
-
-       const std::string& configFileName   () const { return _configFileName;      }  ///< returns path to config file
-       const std::string& eventDataFileName() const { return _eventDataFileName;   }  ///< returns path to output file
-
-       void setConfigFileName   (const std::string& configFileName   ) { _configFileName    = configFileName;    }  ///< sets path to config file
-       void setEventDataFileName(const std::string& eventDataFileName) { _eventDataFileName = eventDataFileName; }  ///< sets path to output file
-       
-       void boostEvent(upcEvent &e); ///< Boost event from beam CMS to lab system
-
-private:
-
-       std::string      _configFileName;     ///< path to configuration file
-       std::string      _eventDataFileName;  ///< path to output file
-
-       starlight*       _starlight;         ///< pointer to starlight instance
-       inputParameters* _inputParameters;   ///< pointer to parameter instance
-       unsigned int     _nmbEventsTot;      ///< total number of events to generate (taken from configuration file)
-       unsigned int     _nmbEventsPerFile;  ///< maximum number of events written to a single file (not yet implemented)
-
-};
-
-
-#endif  // STARLIGHTSTANDALONE_H
diff --git a/STARLIGHT/starlight/src/beam.cpp~ b/STARLIGHT/starlight/src/beam.cpp~
deleted file mode 100644 (file)
index 286825e..0000000
+++ /dev/null
@@ -1,86 +0,0 @@
-///////////////////////////////////////////////////////////////////////////
-//
-//    Copyright 2010
-//
-//    This file is part of starlight.
-//
-//    starlight 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 3 of the License, or
-//    (at your option) any later version.
-//
-//    starlight 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 starlight. If not, see <http://www.gnu.org/licenses/>.
-//
-///////////////////////////////////////////////////////////////////////////
-//
-// File and Version Information:
-// $Rev:: 164                         $: revision of last commit
-// $Author:: odjuvsla                 $: author of last commit
-// $Date:: 2013-10-06 16:18:08 +0200 #$: date of last commit
-//
-// Description:
-//
-//
-//
-///////////////////////////////////////////////////////////////////////////
-
-
-#include <iostream>
-#include <fstream>
-#include <cmath>
-
-#include "starlightconstants.h"
-#include "inputParameters.h"
-#include "reportingUtils.h"
-#include "bessel.h"
-#include "beam.h"
-
-
-using namespace std;
-using namespace starlightConstants;
-
-//______________________________________________________________________________
-beam::beam(const int              Z,
-           const int              A,
-           const double           bdeuteron,
-           const bool             dAuCoherentProduction)
-       : nucleus(Z,
-                 A,
-                 bdeuteron,
-                 dAuCoherentProduction)
-       ,_beamLorentzGamma(inputParametersInstance.beamLorentzGamma())
-{
-}
-
-
-//______________________________________________________________________________
-beam::~beam()
-{ }
-
-
-//______________________________________________________________________________
-double beam::photonFlux(const double impactparameter, 
-                        const double photonEnergy) const
-{
-  // function for the calculation of the "photon density".
-  // photonFlux = number of photons / (energy * area)
-  // assume beta = 1 and gamma >> 1, i.e. neglect the (1 / gamma^2) * K_0(x) term
-  
-  const double X
-         = (impactparameter * photonEnergy) / (_beamLorentzGamma * hbarc);
-  if (X <= 0) 
-         printWarn << "X = " << X << endl;
-  
-  const double factor1 = (double(Z() * Z()) * alpha) / (pi * pi);  
-  const double factor2 = 1. / (photonEnergy * impactparameter * impactparameter);
-  const double bessel  = bessel::dbesk1(X);
-  const double factor3 = X * X * bessel * bessel;
-
-  return factor1 * factor2 * factor3;
-}
diff --git a/STARLIGHT/starlight/src/beambeamsystem.cpp~ b/STARLIGHT/starlight/src/beambeamsystem.cpp~
deleted file mode 100644 (file)
index 41d0fd9..0000000
+++ /dev/null
@@ -1,669 +0,0 @@
-///////////////////////////////////////////////////////////////////////////
-//
-//    Copyright 2010
-//
-//    This file is part of starlight.
-//
-//    starlight 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 3 of the License, or
-//    (at your option) any later version.
-//
-//    starlight 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 starlight. If not, see <http://www.gnu.org/licenses/>.
-//
-///////////////////////////////////////////////////////////////////////////
-//
-// File and Version Information:
-// $Rev:: 164                         $: revision of last commit
-// $Author:: odjuvsla                 $: author of last commit
-// $Date:: 2013-10-06 16:18:08 +0200 #$: date of last commit
-//
-// Description:
-//
-//
-//
-///////////////////////////////////////////////////////////////////////////
-
-
-#include <iostream>
-#include <fstream>
-#include <cmath>
-
-#include "inputParameters.h"
-#include "reportingUtils.h"
-#include "starlightconstants.h"
-#include "bessel.h"
-#include "beambeamsystem.h"
-
-
-using namespace std;
-using namespace starlightConstants;
-
-
-//______________________________________________________________________________
-beamBeamSystem::beamBeamSystem(const beam&            beam1,
-                               const beam&            beam2)
-  : _beamLorentzGamma(inputParametersInstance.beamLorentzGamma()),
-    _beamBreakupMode (inputParametersInstance.beamBreakupMode()),
-    _beam1           (beam1),
-    _beam2           (beam2),
-    _breakupProbabilities(0),
-    _breakupImpactParameterStep(1.007),
-    _breakupCutOff(10e-6)
-{ 
-  init();
-}
-  
-
-
-
-//______________________________________________________________________________
-beamBeamSystem::beamBeamSystem()
-       : _beamLorentzGamma(inputParametersInstance.beamLorentzGamma()),
-         _beamBreakupMode (inputParametersInstance.beamBreakupMode()),
-         _beam1           (inputParametersInstance.beam1Z(),
-                           inputParametersInstance.beam1A(),
-                           inputParametersInstance.deuteronSlopePar(),
-                           inputParametersInstance.coherentProduction()),
-         _beam2           (inputParametersInstance.beam2Z(),
-                           inputParametersInstance.beam2A(),
-                           inputParametersInstance.deuteronSlopePar(),
-                           inputParametersInstance.coherentProduction()),
-         _breakupProbabilities(0),
-         _breakupImpactParameterStep(1.007),
-         _breakupCutOff(10e-10)
-{
-  init();
-}
-
-
-
-//______________________________________________________________________________
-beamBeamSystem::~beamBeamSystem()
-{ }
-
-void beamBeamSystem::init()
-{
-   // Calculate beam gamma in CMS frame
-   double rap1 = acosh(inputParametersInstance.beam1LorentzGamma());
-   double rap2 = -acosh(inputParametersInstance.beam2LorentzGamma());
-   
-   _cmsBoost = (rap1+rap2)/2.;
-
-   _beamLorentzGamma = cosh((rap1-rap2)/2);
-   _beam1.setBeamLorentzGamma(_beamLorentzGamma);
-   _beam2.setBeamLorentzGamma(_beamLorentzGamma);
-   
-   generateBreakupProbabilities();
-}
-//______________________________________________________________________________
-double
-beamBeamSystem::probabilityOfBreakup(const double D) const
-{
-       
-       double bMin = (_beam1.nuclearRadius()+_beam2.nuclearRadius())/2.;
-       double pOfB = 0.; // PofB = 1 means that there will be a UPC event and PofB = 0 means no UPC
-
-       // Do pp here
-       if ((_beam1.Z() == 1) && (_beam1.A() == 1) && (_beam2.Z() == 1) && (_beam2.A() == 1)) {  
-               double ppslope=19.8;
-               double GammaProfile = exp(-D * D / (2. * hbarc * hbarc * ppslope));
-               pOfB = (1. - GammaProfile) * (1. - GammaProfile);
-               // if (D < 2. * _beam1.nuclearRadius())
-               //      //Should be the total of RNuc1+Rnuc2,used only beam #1
-               //      PofB = 0.0;
-               // else
-               //      PofB = 1.0;
-               return pOfB;
-       }
-       else if ( ( (_beam1.A() == 1) && (_beam2.A() != 1) ) || ((_beam1.A() != 1) && (_beam2.A() == 1)) ) {  
-         // This is pA
-          if( _beam1.A() == 1 ){ 
-            bMin = _beam2.nuclearRadius() + 0.7; 
-          }else if( _beam2.A() == 1 ){ 
-            bMin = _beam1.nuclearRadius() + 0.7; 
-          }else{
-            cout<<"Some logical problem here!"<<endl;
-          }
-          if( D > bMin )pOfB=1.0; 
-          return pOfB;
-       }
-
-       //use the lookup table and return...
-       pOfB = 1.;
-       if (D > 0.0) {             
-               //Now we must determine which step number in d corresponds to this D,
-               // and use appropiate Ptot(D_i)
-               //int i = (int)(log(D / Bmin) / log(1.01));
-               int i = (int)(log(D / bMin) / log(_breakupImpactParameterStep));
-               if (i <= 0)
-                       pOfB = _breakupProbabilities[0];
-               else{
-                       if (i >= int(_breakupProbabilities.size()-1))
-                               pOfB = _breakupProbabilities[_breakupProbabilities.size()-1];
-                       else {
-                               // const double DLow = Bmin * pow((1.01), i);
-                               const double DLow = bMin * pow((_breakupImpactParameterStep), i);
-                               // const double DeltaD = 0.01 * DLow;
-                               const double DeltaD = (_breakupImpactParameterStep-1) * DLow;
-                               const double DeltaP = _breakupProbabilities[i + 1] - _breakupProbabilities[i];
-                               pOfB   = _breakupProbabilities[i] + DeltaP * (D - DLow) / DeltaD;
-                       }
-               }
-       }
-
-       return pOfB;
-}
-
-void
-beamBeamSystem::generateBreakupProbabilities()
-{
-    // Step = 1.007;//.01; //We will multiplicateively increase Biter by 1%
-    
-    
-    double bMin = (_beam1.nuclearRadius()+_beam2.nuclearRadius())/2.;
-    
-    
-    if ((_beam1.Z() != 1) && (_beam1.A() != 1) && (_beam2.Z() != 1) && _beam2.A() != 1) {
-
-        if (_beamBreakupMode == 1)
-            printInfo << "Hard Sphere Break criteria. b > " << 2. * _beam1.nuclearRadius() << endl;
-        if (_beamBreakupMode == 2)
-            printInfo << "Requiring XnXn [Coulomb] breakup. " << endl;
-        if (_beamBreakupMode == 3)
-            printInfo << "Requiring 1n1n [Coulomb only] breakup. " << endl;
-        if (_beamBreakupMode == 4)
-            printInfo << "Requiring both nuclei to remain intact. " << endl;
-        if (_beamBreakupMode == 5)
-            printInfo << "Requiring no hadronic interactions. " << endl;
-        if (_beamBreakupMode == 6)
-            printInfo << "Requiring breakup of one or both nuclei. " << endl;
-        if (_beamBreakupMode == 7)
-            printInfo << "Requiring breakup of one nucleus (Xn,0n). " << endl;
-
-        //pp may cause segmentation fault in here and it does not use this...
-       double pOfB = 0;
-       double b = bMin;
-        double totRad = _beam1.nuclearRadius()+_beam2.nuclearRadius();
-        
-       while(1)
-       {
-            
-            if(_beamBreakupMode != 5)
-            {
-                if(b > (totRad*1.5))
-                {
-                    if(pOfB<_breakupCutOff)
-                    {
-//                         std::cout << "Break off b: " << b << std::endl;
-//                         std::cout << "Number of PofB bins: " << _breakupProbabilities.size() << std::endl;
-                        break;
-                    }
-                }
-            }
-            else
-            {
-                if((1-pOfB)<_breakupCutOff)
-                {
- //                         std::cout << "Break off b: " << b << std::endl;
-//                         std::cout << "Number of PofB bins: " << _breakupProbabilities.size() << std::endl;
-                        break;
-                }
-            }
-//             std::cout << 1-pOfBreakup << std::endl;
-//             _pHadronBreakup = 0;
-//             _pPhotonBreakup = 0;
-
-//             double pHadronBreakup = probabilityOfHadronBreakup(b);
-            probabilityOfHadronBreakup(b);
-            //moved gammatarg into photonbreakup
-//             double pPhotonBreakup = probabilityOfPhotonBreakup(b, _beamBreakupMode);
-            probabilityOfPhotonBreakup(b, _beamBreakupMode);
-
-            //What was probability of photonbreakup depending upon mode selection,
-            // is now done in the photonbreakupfunction
-            if (_beamBreakupMode == 1) {
-                if (b >_beam1.nuclearRadius()+_beam2.nuclearRadius())  // symmetry
-                    _pHadronBreakup = 0;
-                else
-                    _pHadronBreakup = 999.;
-            }
-            
-            b *= _breakupImpactParameterStep;
-           pOfB = exp(-1 * _pHadronBreakup) * _pPhotonBreakup;
-            _breakupProbabilities.push_back(pOfB);
-        }
-    }
-    else if (((_beam1.Z() == 1) && (_beam1.A() == 1)) || ((_beam2.Z() == 1) && (_beam2.A() == 1))) {  
-      
-      double pOfB = 0;
-      double b = bMin;
-      double totRad = _beam1.nuclearRadius()+_beam2.nuclearRadius();
-      
-      while(1)
-       {
-            if(_beamBreakupMode != 5)
-            {
-                if(b > (totRad*1.5))
-                {
-                    if(pOfB<_breakupCutOff)
-                    {
-//                         std::cout << "Break off b: " << b << std::endl;
-                        break;
-                    }
-                }
-            }
-            else
-            {
-                if((1-pOfB)<_breakupCutOff)
-                {
-//                         std::cout << "Break off b: " << b << std::endl;
-                        break;
-                }
-            }
-         _beam1.Z() > 1 ? pOfB = exp(-7.35*_beam1.thickness(b)) :
-                          pOfB = exp(-7.35*_beam2.thickness(b));
-         _breakupProbabilities.push_back(pOfB);
-            b *= _breakupImpactParameterStep;
-        }
-    }
-
-    
-}
-
-//______________________________________________________________________________
-double
-beamBeamSystem::probabilityOfHadronBreakup(const double impactparameter)
-{
-       //      double pbreakup =0.; 
-       //probability of hadron breakup, 
-       //this is what is returned when the function is called
-       double gamma = _beamLorentzGamma; 
-       //input for gamma_em
-       //this will need to be StarlightInputParameters::gamma or whatever
-       double b = impactparameter;
-       int ap = _beam1.A();  
-       //Notice this is  taking from nucleus 1.Still assuming symmetric system?
-
-       static int IFIRSTH = 0;
-       static double DELL=0., DELR=0., SIGNN=0., R1=0., A1=0., R2=0., RHO1=0.;
-       static double RHO2=0., NZ1=0., NZ2=0., NR1=0., NR2=0.,RR1=0., NY=0., NX=0.;
-       static double AN1=0., AN2=0.;
-       double delo=0.,RSQ=0.,Z1=0.,Y=0.,X=0.,XB=0.,RPU=0.,IRUP=0.,RTU=0.;
-       double IRUT=0.,T1=0.,T2=0.;
-       static double DEN1[20002], DEN2[20002];
-       if (IFIRSTH != 0) goto L100;
-       //Initialize
-       //Integration delta x, delta z
-       IFIRSTH = 1;
-       DELL   = .05;
-       DELR   = .01;
-
-       //use two sigma_NN's. 52mb at rhic 100gev/beam, 88mb at LHC 2.9tev/beam, gamma is in cm system
-       SIGNN = 5.2;
-       if ( gamma > 500. ) SIGNN = 8.8;
-       //use parameter from Constants
-       R1 = ( _beam1.nuclearRadius());  //remember _beam2? better way to do this generically
-       A1 = 0.535; //This is woodsaxonskindepth?
-       //write(6,12)r1,a1,signn  Here is where we could probably set this up asymmetrically R2=_beam2.nuclearRadius() and RHO2=ap2=_beam2.A()
-       R2 = R1;
-       RHO1 = ap;
-       RHO2 = RHO1;
-       NZ1  = ((R1+5.)/DELR);
-       NR1  = NZ1;
-       NZ2  = ((R2+5.)/DELR);
-       NR2  = NZ2;
-       RR1  = -DELR;
-       NY   = ((R1+5.)/DELL);
-       NX   = 2*NY;
-       for ( int IR1 = 1; IR1 <= NR1; IR1++) {
-               DEN1[IR1] = 0.;
-               RR1       = RR1+DELR;
-               Z1        = -DELR/2;
-
-               for ( int IZ1 = 1; IZ1 <= NZ1; IZ1++) {
-                       Z1  = Z1+DELR;
-                       RSQ = RR1*RR1+Z1*Z1;
-                       DEN1[IR1] = DEN1[IR1]+1./(1.+exp((sqrt(RSQ)-R1)/A1));
-               }//for(IZ)
-
-               DEN1[IR1] = DEN1[IR1]*2.*DELR;
-               DEN2[IR1] = DEN1[IR1];
-       }//for(IR)
-
-       AN1 = 0.;
-       RR1 = 0.;
-    
-       for ( int IR1 =1; IR1 <= NR1; IR1++) {
-               RR1 = RR1+DELR;
-               AN1 = AN1+RR1*DEN1[IR1]*DELR*2.*3.141592654;
-       }//for(IR)
-        
-       AN2 = AN1; //This will also probably need to be changed?
-
-       delo = .05;
-       //.1 to turn mb into fm^2
-       //Calculate breakup probability here
- L100:
-       _pHadronBreakup = 0.;
-       if ( b > 25. ) return _pHadronBreakup;
-       Y = -.5*DELL;
-       for ( int IY = 1; IY <= NY; IY++) {
-               Y = Y+DELL;
-               X = -DELL*float(NY+1);
-
-               for ( int IX = 1; IX <=NX; IX++) {
-                       X = X+DELL;
-                       XB = b-X;
-                       RPU = sqrt(X*X+Y*Y);
-                       IRUP = (RPU/DELR)+1;
-                       RTU  = sqrt(XB*XB+Y*Y);
-                       IRUT = (RTU/DELR)+1;
-                       T1   = DEN2[(int)IRUT]*RHO2/AN2;
-                       T2   = DEN1[(int)IRUP]*RHO1/AN1;
-                       //Eq.6 BCW, Baltz, Chasman, White, Nucl. Inst. & Methods A 417, 1 (1998)
-                       _pHadronBreakup=_pHadronBreakup+2.*T1*(1.-exp(-SIGNN*T2))*DELL*DELL;
-               }//for(IX)
-       }//for(IY)
-
-       return _pHadronBreakup;
-}
-
-
-//______________________________________________________________________________
-double
-beamBeamSystem::probabilityOfPhotonBreakup(const double impactparameter, const int mode)
-{
-       static double ee[10001], eee[162], se[10001];
-
-       //double gamma_em=108.4;  //This will be an input value.
-       _pPhotonBreakup =0.;   //Might default the probability with a different value?
-       double b = impactparameter;
-       int zp = _beam1.Z();  //What about _beam2? Generic approach?
-       int ap = _beam1.A();
-       
-       //Was initialized at the start of the function originally, been moved inward.
-       double pxn=0.;
-       double p1n=0.;
-
-       //Used to be done prior to entering the function. Done properly for assymetric?
-       double gammatarg = 2.*_beamLorentzGamma*_beamLorentzGamma-1.;   
-       double omaxx =0.;
-       //This was done prior entering the function as well
-       if (_beamLorentzGamma > 500.){
-               omaxx=1.E10;
-       }
-       else{
-               omaxx=1.E7;
-       }
-
-
-       double e1[23]= {0.,103.,106.,112.,119.,127.,132.,145.,171.,199.,230.,235.,
-                       254.,280.,300.,320.,330.,333.,373.,390.,420.,426.,440.};
-       double s1[23]= {0.,12.0,11.5,12.0,12.0,12.0,15.0,17.0,28.0,33.0,
-                       52.0,60.0,70.0,76.0,85.0,86.0,89.0,89.0,75.0,76.0,69.0,59.0,61.0};
-       double e2[12]={0.,2000.,3270.,4100.,4810.,6210.,6600.,
-                      7790.,8400.,9510.,13600.,16400.};
-       double s2[12]={0.,.1266,.1080,.0805,.1017,.0942,.0844,.0841,.0755,.0827,
-                      .0626,.0740};
-       double e3[29]={0.,26.,28.,30.,32.,34.,36.,38.,40.,44.,46.,48.,50.,52.,55.,
-                      57.,62.,64.,66.,69.,72.,74.,76.,79.,82.,86.,92.,98.,103.};
-       double s3[29]={0.,30.,21.5,22.5,18.5,17.5,15.,14.5,19.,17.5,16.,14.,
-                      20.,16.5,17.5,17.,15.5,18.,15.5,15.5,15.,13.5,18.,14.5,15.5,12.5,13.,
-                      13.,12.};
-       static double sa[161]={0.,0.,.004,.008,.013,.017,.021,.025,.029,.034,.038,.042,.046,
-                              .051,.055,.059,.063,.067,.072,.076,.08,.085,.09,.095,.1,.108,.116,
-                              .124,.132,.14,.152,.164,.176,.188,.2,.22,.24,.26,.28,.3,.32,.34,
-                              .36,.38,.4,.417,.433,.450,.467,.483,.5,.51,.516,.52,.523,.5245,
-                              .525,.5242,
-                              .5214,.518,.512,.505,.495,.482,.469,.456,.442,.428,.414,.4,.386,
-                              .370,.355,.34,.325,.310,.295,.280,.265,.25,.236,.222,.208,.194,
-                              .180,.166,
-                              .152,.138,.124,.11,.101,.095,.09,.085,.08,.076,.072,.069,.066,
-                              .063,.06,.0575,.055,.0525,.05,.04875,.0475,.04625,.045,.04375,
-                              .0425,.04125,.04,.03875,.0375,.03625,.035,.03375,.0325,.03125,.03,
-                              .02925,.0285,.02775,.027,.02625,.0255,.02475,.024,.02325,.0225,
-                              .02175,.021,.02025,.0195,.01875,.018,.01725,.0165,.01575,.015,
-                              .01425,.0135,.01275,.012,.01125,.0105,.00975,.009,.00825,.0075,
-                              .00675,.006,.00525,.0045,.00375,.003,.00225,.0015,.00075,0.};
-
-
-
-       double sen[161]={0.,0.,.012,.025,.038,.028,.028,.038,.035,.029,.039,.035,
-                        .038,.032,.038,.041,.041,.049,.055,.061,.072,.076,.070,.067,
-                        .080,.103,.125,.138,.118,.103,.129,.155,.170,.180,.190,.200,
-                        .215,.250,.302,.310,.301,.315,.330,.355,.380,.400,.410,.420,
-                        .438,.456,.474,.492,.510,.533,.556,.578,.6,.62,.63,.638,
-                        .640,.640,.637,.631,.625,.618,.610,.600,.580,.555,.530,.505,
-                        .480,.455,.435,.410,.385,.360,.340,.320,.300,.285,.270,.255,
-                        .240,.225,.210,.180,.165,.150,.140,.132,.124,.116,.108,.100,
-                        .092,.084,.077,.071,.066,.060,.055,.051,.048,.046,.044,.042,
-                        .040,.038,.036,.034,.032,.030,.028,.027,.026,.025,.025,.025,
-                        .024,.024,.024,.024,.024,.023,.023,.023,.023,.023,.022,.022,
-                        .022,.022,.022,.021,.021,.021,.020,.020,
-                        .020,.019,.018,.017,.016,.015,.014,.013,.012,.011,.010,.009,
-                        .008,.007,.006,.005,.004,.003,.002,.001,0.};
-
-       // gammay,p gamma,n of Armstrong begin at 265 incr 25
-
-
-       double sigt[160]={0.,.4245,.4870,.5269,.4778,.4066,.3341,.2444,.2245,.2005,
-                         .1783,.1769,.1869,.1940,.2117,.2226,.2327,.2395,.2646,.2790,.2756,
-                         .2607,.2447,.2211,.2063,.2137,.2088,.2017,.2050,.2015,.2121,.2175,
-                         .2152,.1917,.1911,.1747,.1650,.1587,.1622,.1496,.1486,.1438,.1556,
-                         .1468,.1536,.1544,.1536,.1468,.1535,.1442,.1515,.1559,.1541,.1461,
-                         .1388,.1565,.1502,.1503,.1454,.1389,.1445,.1425,.1415,.1424,.1432,
-                         .1486,.1539,.1354,.1480,.1443,.1435,.1491,.1435,.1380,.1317,.1445,
-                         .1375,.1449,.1359,.1383,.1390,.1361,.1286,.1359,.1395,.1327,.1387,
-                         .1431,.1403,.1404,.1389,.1410,.1304,.1363,.1241,.1284,.1299,.1325,
-                         .1343,.1387,.1328,.1444,.1334,.1362,.1302,.1338,.1339,.1304,.1314,
-                         .1287,.1404,.1383,.1292,.1436,.1280,.1326,.1321,.1268,.1278,.1243,
-                         .1239,.1271,.1213,.1338,.1287,.1343,.1231,.1317,.1214,.1370,.1232,
-                         .1301,.1348,.1294,.1278,.1227,.1218,.1198,.1193,.1342,.1323,.1248,
-                         .1220,.1139,.1271,.1224,.1347,.1249,.1163,.1362,.1236,.1462,.1356,
-                         .1198,.1419,.1324,.1288,.1336,.1335,.1266};
-
-
-       double sigtn[160]={0.,.3125,.3930,.4401,.4582,.3774,.3329,.2996,.2715,.2165,
-                          .2297,.1861,.1551,.2020,.2073,.2064,.2193,.2275,.2384,.2150,.2494,
-                          .2133,.2023,.1969,.1797,.1693,.1642,.1463,.1280,.1555,.1489,.1435,
-                          .1398,.1573,.1479,.1493,.1417,.1403,.1258,.1354,.1394,.1420,.1364,
-                          .1325,.1455,.1326,.1397,.1286,.1260,.1314,.1378,.1353,.1264,.1471,
-                          .1650,.1311,.1261,.1348,.1277,.1518,.1297,.1452,.1453,.1598,.1323,
-                          .1234,.1212,.1333,.1434,.1380,.1330,.12,.12,.12,.12,.12,.12,.12,.12,
-                          .12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,
-                          .12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,
-                          .12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,
-                          .12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,
-                          .12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12,.12};
-       //89*.12};
-
-
-
-       static int IFIRSTP=0;
-
-
-       //  Initialization needed?
-       //double hbar=197.3;
-       //double pi=3.141592654;
-
-       // added
-       double si1=0, g1 =0,   o1=0;
-       int   ne = 0, ij =0;
-       double delo=0, omax =0, gk1m=0;
-       static double scon=0., zcon=0.,o0=0.;
-
-
-       double x=0,y=0,eps=0,eta=0,em=0,exx=0,s=0,ictr=0,pom=0,vec=0,gk1=0;
-
-       //  maximum energy for GDR dissocation (in target frame, in MeV)
-
-       double omax1n=24.01;
-
-       if (IFIRSTP != 0) goto L100;
-
-       IFIRSTP=1;
-
-
-       //This whole thing is dependenant on gold or lead....Might need to expand
-       if (zp == 79)
-               {
-
-
-                       ap=197;
-                       si1=540.;
-                       g1=4.75;
-
-                       // peak and minimum energies for GDR excitation (in MeV)
-                       o1=13.70;
-                       o0=8.1;
-               }
-       else
-               {
-                       zp=82;   //assumed to be lead
-                       ap=208;
-                       si1=640.;
-                       g1=4.05;
-                       o1=13.42;
-                       o0=7.4;
-                       for(int j=1;j<=160;j++)
-                               {
-
-                                       sa[j]=sen[j];
-                               }
-               }
-       //Part II of initialization
-       delo = .05;
-       //.1 to turn mb into fm^2
-       scon = .1*g1*g1*si1;
-       zcon = zp/(gammatarg*( pi)*( 
-                                                       hbarcmev))*zp/(gammatarg*( pi)*
-                            ( hbarcmev))/137.04;
-
-       //single neutron from GDR, Veyssiere et al. Nucl. Phys. A159, 561 (1970)
-       for ( int i = 1; i <= 160; i++) {
-               eee[i] = o0+.1*(i-1);
-               sa[i]  = 100.*sa[i];
-       }
-       //See Baltz, Rhoades-Brown, and Weneser, Phys. Rev. E 54, 4233 (1996) 
-       //for details of the following photo cross-sections
-       eee[161]=24.1;
-       ne=int((25.-o0)/delo)+1;
-       //GDR any number of neutrons, Veyssiere et al., Nucl. Phys. A159, 561 (1970)
-       for ( int i = 1; i <= ne; i++ ) {
-               ee[i] = o0+(i-1)*delo;
-               //cout<<" ee 1 "<<ee[i]<<"  "<<i<<endl;
-
-               se[i] = scon*ee[i]*ee[i]/(((o1*o1-ee[i]*ee[i])*(o1*o1-ee[i]*ee[i]))
-                                         +ee[i]*ee[i]*g1*g1);
-       }
-       ij = ne;   //Risky?
-       //25-103 MeV, Lepretre, et al., Nucl. Phys. A367, 237 (1981)
-       for ( int j = 1; j <= 27; j++ ) {
-               ij = ij+1;
-               ee[ij] = e3[j];
-               //cout<<" ee 2 "<<ee[ij]<<"  "<<ij<<endl;
-
-               se[ij] = .1*ap*s3[j]/208.;
-       }
-       //103-440 MeV, Carlos, et al., Nucl. Phys. A431, 573 (1984)
-       for ( int j = 1; j <= 22; j++ ) {
-               ij = ij+1;
-               ee[ij] = e1[j];
-               //cout<<" ee 3 "<<ee[ij]<<"  "<<ij<<endl;
-               se[ij] = .1*ap*s1[j]/208.;
-       }
-       //440 MeV-2 GeV Armstrong et al.
-       for ( int j = 9; j <= 70; j++) {
-               ij = ij+1;
-               ee[ij] = ee[ij-1]+25.;
-               //cout<<" ee 4 "<<ee[ij]<<"  "<<ij<<endl;
-               se[ij] = .1*(zp*sigt[j]+(ap-zp)*sigtn[j]);
-       }
-       //2-16.4 GeV Michalowski; Caldwell
-       for ( int j = 1; j <= 11; j++) {
-               ij = ij+1;
-               ee[ij] = e2[j];
-               //cout<<" ee 5 "<<ee[ij]<<"   "<<ij<<endl;
-               se[ij] = .1*ap*s2[j];
-       }
-       //Regge paramteres
-       x = .0677;
-       y = .129;
-       eps = .0808;
-       eta = .4525;
-       em = .94;
-       exx = pow(10,.05);
-
-       //Regge model for high energy
-       s = .002*em*ee[ij];
-       //make sure we reach LHC energies
-       ictr = 100;
-       if ( gammatarg > (2.*150.*150.)) ictr = 150;
-       for ( int j = 1; j <= ictr; j++ ) {
-               ij = ij+1;
-               s = s*exx;
-               ee[ij] = 1000.*.5*(s-em*em)/em;
-               //cout<<" ee 6 "<<ee[ij]<<"   "<<ij<<endl;
-               pom = x*pow(s,eps);
-               vec = y*pow(s,(-eta));
-               se[ij] = .1*.65*ap*(pom+vec);
-       }
-       ee[ij+1] = 99999999999.;
-       //done with initaliation
-       //write(6,99)o0;
-       //clear counters for 1N, XN
- L100:
-
-       p1n = 0.;
-       pxn = 0.;
-       //start XN calculation
-       //what's the b-dependent highest energy of interest?
-
-       omax = min(omaxx,4.*gammatarg*( hbarcmev)/b);
-       if ( omax < o0 ) return _pPhotonBreakup;
-       gk1m = bessel::dbesk1(ee[1]*b/(( hbarcmev)*gammatarg));
-       int k = 2;
- L212:
-       if (ee[k] < omax ) {
-               gk1 = bessel::dbesk1(ee[k]*b/(( hbarcmev)*gammatarg));
-               //Eq. 3 of BCW--NIM in Physics Research A 417 (1998) pp1-8:
-               pxn=pxn+zcon*(ee[k]-ee[k-1])*.5*(se[k-1]*ee[k-1]*gk1m*gk1m+se[k]*ee[k]*gk1*gk1);
-               k = k + 1;
-               gk1m = gk1;
-               goto L212;
-       }
-       //one neutron dissociation
-       omax = min(omax1n,4.*gammatarg*( hbarcmev)/b);
-       gk1m = bessel::dbesk1(eee[1]*b/(( hbarcmev)*gammatarg));
-       k = 2;
- L102:
-       if (eee[k] < omax ) {
-               gk1 = bessel::dbesk1(eee[k]*b/(( hbarcmev)*gammatarg));
-               //Like Eq3 but with only the one neutron out GDR photo cross section input
-               p1n = p1n+zcon*(eee[k]-eee[k-1])*.5*(sa[k-1]*eee[k-1]*gk1m*gk1m+sa[k]*eee[k]*gk1*gk1);
-               k = k+1;
-               gk1m = gk1;
-               goto L102;
-       }
-
-
-       //This used to be done externally, now it is done internally.
-       if (( mode) == 1) _pPhotonBreakup = 1.;
-       if (( mode) == 2) _pPhotonBreakup = (1-exp(-1*pxn))*(1-exp(-1*pxn));
-       if (( mode) == 3) _pPhotonBreakup = (p1n*exp(-1*pxn))*(p1n*exp(-1*pxn));
-       if (( mode) == 4) _pPhotonBreakup = exp(-2*pxn);
-       if (( mode) == 5) _pPhotonBreakup = 1.;
-       if (( mode) == 6) _pPhotonBreakup = (1. - exp(-2.*pxn));
-       if (( mode) == 7) _pPhotonBreakup = 2.*exp(-pxn)*(1.-exp(-pxn));
-
-       //cout<<pxn<<" "<<zcon<<" "<<ee[k]<<" "<<se[k-1]<<" "<<gk1m<<"  "<<gk1<<"  "<<k<<"  "<<ee[k+1]<< "  "<<b<< endl;
-
-       return _pPhotonBreakup;
-}
diff --git a/STARLIGHT/starlight/src/gammagammaleptonpair.cpp~ b/STARLIGHT/starlight/src/gammagammaleptonpair.cpp~
deleted file mode 100644 (file)
index eb93344..0000000
+++ /dev/null
@@ -1,848 +0,0 @@
-///////////////////////////////////////////////////////////////////////////
-//
-//    Copyright 2010
-//
-//    This file is part of starlight.
-//
-//    starlight 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 3 of the License, or
-//    (at your option) any later version.
-//
-//    starlight 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 starlight. If not, see <http://www.gnu.org/licenses/>.
-//
-///////////////////////////////////////////////////////////////////////////
-//
-// File and Version Information:
-// $Rev:: 164                         $: revision of last commit
-// $Author:: odjuvsla                 $: author of last commit
-// $Date:: 2013-10-06 16:18:08 +0200 #$: date of last commit
-//
-// Description:
-//    Nystrand 220710
-//    Fixed bug which gave incorrect minv distribution in gammagammaleptonpair.
-//    Moved loop over W and Y from pickw to twoLeptonCrossSection in
-//    gammagammaleptonpair to speed up event generation.
-//    Changed to Woods-Saxon radius in twophotonluminosity to be consistent
-//    with old starligt.
-//
-//
-///////////////////////////////////////////////////////////////////////////
-
-
-#include <iostream>
-#include <fstream>
-#include <cmath>
-#include <vector>
-
-#include "starlightconstants.h"
-#include "gammagammaleptonpair.h"
-
-
-using namespace std;
-
-
-//_____________________________________________________________________________
-Gammagammaleptonpair::Gammagammaleptonpair(beamBeamSystem& bbsystem)
-: eventChannel(bbsystem)
-{
-    //Initialize randomgenerator with our seed.
-    cout<<"Randy in leptonpair construction: "<<randyInstance.Rndom()<<endl;
-    //Storing inputparameters into protected members for use
-    _GGlepInputnumw=inputParametersInstance.nmbWBins();
-    _GGlepInputnumy=inputParametersInstance.nmbRapidityBins();
-    _GGlepInputpidtest=inputParametersInstance.prodParticleType();
-    _GGlepInputGamma_em=inputParametersInstance.beamLorentzGamma();
-    //Let us read in the luminosity tables
-    read();
-    //Now we will calculate the crosssection
-    twoLeptonCrossSection();
-    //If it is a tauon, calculate its tables
-    if(inputParametersInstance.prodParticleId()==starlightConstants::TAUON) calculateTable();
-}
-
-
-//______________________________________________________________________________
-Gammagammaleptonpair::~Gammagammaleptonpair()
-{ }
-
-
-//______________________________________________________________________________
-void Gammagammaleptonpair::twoLeptonCrossSection()
-{
-    //This function calculates section for 2-particle decay. For reference, see STAR Note 243, Eq. 9.
-    //calculate the 2-lepton differential cross section
-    //the 100 is to convert to barns
-    //the 2 is to account for the fact that we integrate only over one half of the rapidity range
-    //Multiply all _Farray[] by _f_max
-
-    for(int i=0;i<_GGlepInputnumw;i++)
-    {
-       for(int j=0;j<_GGlepInputnumy;j++)
-       {
-           // _sigmax[i][j]=2.*Gammagammaleptonpair::twoMuonCrossSection(_Warray[i])*_f_max*_Farray[i][j]/100.;
-            _sigmax[i][j]=twoMuonCrossSection(_Warray[i])*_f_max*_Farray[i][j]/(100.*_Warray[i]);
-       }
-    }
-    //calculate the total two-lepton cross section
-    double sigmasum =0.;
-    for(int i=0;i<_GGlepInputnumw-1;i++)
-    {
-       for(int j=0;j<_GGlepInputnumy-1;j++)
-       {
-           // _sigmaSum = _sigmaSum +2.*((_sigmax[i][j]+_sigmax[i+1][j]+_sigmax[i][j+1]+_sigmax[i+1][j+1])/4.*(_Yarray[j+1]-_Yarray[j])*(_Warray[i+1]-_Warray[i])/((_Warray[i+1]+_Warray[i])/2.));
-           // _sigmaSum = _sigmaSum +((_sigmax[i][j]+_sigmax[i+1][j]+_sigmax[i][j+1]+_sigmax[i+1][j+1])/4.*(_Yarray[j+1]-_Yarray[j])*(_Warray[i+1]-_Warray[i])/((_Warray[i+1]+_Warray[i])/2.));
-          sigmasum = sigmasum +(_sigmax[i][j]+_sigmax[i+1][j]+_sigmax[i][j+1]+_sigmax[i+1][j+1])/4.*(_Yarray[j+1]-_Yarray[j])*(_Warray[i+1]-_Warray[i]);
-       }
-    }
-    cout << "The total "<<_GGlepInputpidtest<<" cross-section is: "<<sigmasum<<" barns."<<endl;
-
-    // Do this integration here, once per run rather than once per event (JN 220710) 
-    //integrate sigma down to a function of just w
-
-    double sgf=0.;
-
-    for(int i=0;i<_ReadInputnumw;i++)
-    {
-            _sigofw[i]=0.;
-            for(int j=0;j<_ReadInputnumy-1;j++)
-            {
-                _sigofw[i] = _sigofw[i]+(_Yarray[j+1]-_Yarray[j])*(_sigmax[i][j+1]+_sigmax[i][j])/2.;
-            }
-    }
-
-    //calculate the unnormalized sgfint array
-    _sigfint[0]=0.;
-    for(int i=0;i<_ReadInputnumw-1;i++)
-    {
-        sgf=(_sigofw[i+1]+_sigofw[i])*(_Warray[i+1]-_Warray[i])/2.;
-        _sigfint[i+1]=_sigfint[i]+sgf;
-     }
-
-     //normalize sgfint array
-     _signormw=_sigfint[_ReadInputnumw-1];
-     for(int i=0;i<_ReadInputnumw;i++)
-     {
-          _sigfint[i]=_sigfint[i]/_signormw;
-     }
-    return;
-}
-
-
-//______________________________________________________________________________
-double Gammagammaleptonpair::twoMuonCrossSection(double w)
-{
-    //This function gives the two muon cross section as a function of Y&W. 
-    //Using the formula given in G.Soff et. al Nuclear Equation of State, part B, 579
-    double s=0.,Etest=0.,deltat=0.,xL=0.,sigmuix=0.,alphasquared=0.,hbarcsquared=0.;
-    s = w*w;
-    Etest = 4.*getMass()*getMass()/s;  
-    deltat = s * sqrt(1.-Etest);
-    xL = 2.*log(sqrt(s)/(2.*getMass())+sqrt(1./Etest-1.));
-    alphasquared = starlightConstants::alpha*starlightConstants::alpha;
-    hbarcsquared = starlightConstants::hbarc*starlightConstants::hbarc;
-    sigmuix = 4.*starlightConstants::pi*alphasquared/s*hbarcsquared*((1+Etest-0.5*Etest*Etest)*xL-(1./s+Etest/s)*deltat);
-    if(Etest > 1.) 
-       sigmuix = 0.;
-    return sigmuix;
-}
-
-
-//______________________________________________________________________________
-void Gammagammaleptonpair::pickw(double &w)
-{
-//  This function picks a w for the 2-photon calculation.
-
-    double x=0.,remainarea=0.,remainw=0.,a=0.,b=0.,c=0.;
-    int ivalw=0;
-
-    if(_wdelta != 0)
-    {
-       w=_wdelta;
-       ivalw=_ivalwd;
-       remainw=_remainwd;
-    }
-    else{
-       //deal with the case where sigma is an array
-       //_sigofw is simga integrated over y using a linear interpolation
-       //sigint is the integral of sgfint, normalized
-
-       //pick a random number
-       x = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
-       //compare x and sgfint to find the ivalue which is just less than the random number x
-       for(int i=0;i<_GGlepInputnumw;i++)
-       {
-           if(x > _sigfint[i]) ivalw=i;
-       }
-       //remainder above ivalw
-       remainarea = x - _sigfint[ivalw];
-
-       //figure out what point corresponds to the excess area in remainarea
-       c = -remainarea*_signormw/(_Warray[ivalw+1]-_Warray[ivalw]);
-       b = _sigofw[ivalw];
-       a = (_sigofw[ivalw+1]-_sigofw[ivalw])/2.;
-       if(a==0.)
-       {
-           remainw = -c/b;
-       }
-       else{
-           remainw = (-b+sqrt(b*b-4.*a*c))/(2.*a);
-       }
-       _ivalwd = ivalw;
-       _remainwd = remainw;
-       //calculate the w value
-       w = _Warray[ivalw]+(_Warray[ivalw+1]-_Warray[ivalw])*remainw;
-
-    }
-}
-
-
-//______________________________________________________________________________
-void Gammagammaleptonpair::picky(double &y)
-{
-    // This function picks a y given a W 
-
-    double * sigofy;
-    double * sgfint;
-    sigofy = new double[starlightLimits::MAXYBINS];
-    sgfint = new double[starlightLimits::MAXWBINS];
-       
-    double remainw =0.,remainarea=0.,remainy=0.,a=0.,b=0.,c=0.,sgf=0.,signorm=0.,x=0.;
-    int ivalw=0,ivaly=0;
-
-    ivalw=_ivalwd;
-    remainw=_remainwd;
-    //average over two colums to get y array
-    for(int j=0;j<_GGlepInputnumy;j++)
-    {
-       sigofy[j]=_sigmax[ivalw][j]+(_sigmax[ivalw+1][j]-_sigmax[ivalw][j])*remainw;
-    }
-
-    //calculate the unnormalized sgfint
-    sgfint[0]=0.;
-    for(int j=0;j<_GGlepInputnumy-1;j++)
-    {
-       sgf = (sigofy[j+1]+sigofy[j])/2.;
-       sgfint[j+1]=sgfint[j]+sgf*(_Yarray[j+1]-_Yarray[j]);
-    }
-
-    //normalize the sgfint array
-    signorm = sgfint[_GGlepInputnumy-1];
-
-    for(int j=0;j<_GGlepInputnumy;j++)
-    {
-       sgfint[j]=sgfint[j]/signorm;
-    }
-
-    //pick a random number
-    x = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
-    //compare x and sgfint to find the ivalue which is just less then the random number x
-    for(int i=0;i<_GGlepInputnumy;i++)
-    {
-       if(x > sgfint[i]) ivaly = i;
-    }
-    //remainder above ivaly
-    remainarea = x - sgfint[ivaly];
-
-    //figure what point corresponds to the leftover area in remainarea
-    c = -remainarea*signorm/(_Yarray[ivaly+1]-_Yarray[ivaly]);
-    b = sigofy[ivaly];
-    a = (sigofy[ivaly+1]-sigofy[ivaly])/2.;
-    if(a==0.)
-    {
-       remainy = -c/b;
-    }
-    else{
-       remainy = (-b + sqrt(b*b-4.*a*c))/(2.*a);
-    }
-    //calculate the y value
-    y = _Yarray[ivaly]+(_Yarray[ivaly+1]-_Yarray[ivaly])*remainy;
-    delete[] sigofy;
-    delete[] sgfint;
-}
-
-
-//______________________________________________________________________________
-void Gammagammaleptonpair::pairMomentum(double w,double y,double &E,double &px,double &py,double &pz)
-{
-    //this function calculates px,py,pz,and E given w and y
-
-    double anglepp1=0.,anglepp2=0.,pp1=0.,pp2=0.,E1=0.,E2=0.,signpx=0.,pt=0.;
-
-    //E1 and E2 are for the 2 photons in the CM frame
-    E1 = w*exp(y)/2.;
-    E2 = w*exp(-y)/2.;
-
-    //calculate px and py
-    //to get x and y components-- phi is random between 0 and 2*pi
-    anglepp1 = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
-    anglepp2 = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
-
-    pp1 = pp_1(E1);
-    pp2 = pp_2(E2);
-    px = pp1*cos(2.*starlightConstants::pi*anglepp1)+pp2*cos(2.*starlightConstants::pi*anglepp2);
-    py = pp1*sin(2.*starlightConstants::pi*anglepp1)+pp2*sin(2.*starlightConstants::pi*anglepp2);
-
-    //Compute vector sum Pt=Pt1+Pt2 to find pt for the produced particle
-    pt = sqrt(px*px+py*py);
-
-    //W is the mass of the produced particle (not necessarily on-mass-shell).Now compute its energy and pz
-    E = sqrt(w*w+pt*pt)*cosh(y);
-    pz= sqrt(w*w+pt*pt)*sinh(y);
-    signpx = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
-
-    //pick the z direction
-    //Don't do this anymore since y goes from -ymax to +ymax (JN 15-02-2013)
-    //if(signpx > 0.5) pz = -pz;
-}
-
-
-//______________________________________________________________________________
-double Gammagammaleptonpair::pp_1(double E)
-{
-    // This is for beam 1 
-    // returns on random draw from pp(E) distribution
-    double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.;
-    double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.;
-    int satisfy =0;
-        
-    ereds = (E/_GGlepInputGamma_em)*(E/_GGlepInputGamma_em);
-    //sqrt(3)*E/gamma_em is p_t where the distribution is a maximum
-    Cm = sqrt(3.)*E/_GGlepInputGamma_em;
-    //the amplitude of the p_t spectrum at the maximum
-    singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
-    Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm)));
-        
-    //pick a test value pp, and find the amplitude there
-    x = randyInstance.Rndom();
-    pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius(); 
-    singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
-    test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp)));
-
-    while(satisfy==0){
-       u = randyInstance.Rndom();
-       if(u*Coef <= test)
-       {
-           satisfy =1;
-       }
-       else{
-           x =randyInstance.Rndom();
-           pp = 5*starlightConstants::hbarc/_bbs.beam1().nuclearRadius()*x;
-           singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
-           test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp));
-       }
-    }
-
-    return pp;
-}
-
-//______________________________________________________________________________
-double Gammagammaleptonpair::pp_2(double E)
-{
-
-    // This is for beam 2 
-    //returns on random draw from pp(E) distribution
-    double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.;
-    double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.;
-    int satisfy =0;
-        
-    ereds = (E/_GGlepInputGamma_em)*(E/_GGlepInputGamma_em);
-    //sqrt(3)*E/gamma_em is p_t where the distribution is a maximum
-    Cm = sqrt(3.)*E/_GGlepInputGamma_em;
-    //the amplitude of the p_t spectrum at the maximum
-    singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
-    Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm)));
-        
-    //pick a test value pp, and find the amplitude there
-    x = randyInstance.Rndom(); 
-    pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius(); //Will use nucleus #1 
-    singleformfactorpp1=_bbs.beam2().formFactor(pp*pp+ereds);
-    test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp)));
-
-    while(satisfy==0){
-       u = randyInstance.Rndom(); 
-       if(u*Coef <= test)
-       {
-           satisfy =1;
-       }
-       else{
-           x =randyInstance.Rndom(); 
-           pp = 5*starlightConstants::hbarc/_bbs.beam2().nuclearRadius()*x;
-           singleformfactorpp2=_bbs.beam2().formFactor(pp*pp+ereds); 
-           test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp));
-       }
-    }
-
-    return pp;
-}
-
-
-
-//______________________________________________________________________________
-void Gammagammaleptonpair::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,
-                                    double  ,  // E (unused)
-                                    double  W,
-                                    double  px0, double  py0, double  pz0,
-                                    double& px1, double& py1, double& pz1,
-                                    double& px2, double& py2, double& pz2,
-                                    int&    iFbadevent)
-{
-    //     This routine decays a particle into two particles of mass mdec,
-    //     taking spin into account
-
-    double mdec=0.,E1=0.,E2=0.;
-    double pmag, anglelep[20001];
-    // double ytest=0.,dndtheta;
-    double phi,theta,xtest,Ecm;
-    double betax,betay,betaz;
-    double hirestheta,hirestest,hiresw;  //added from JN->needed precision
-
-    //    set the mass of the daughter particles
-
-    mdec = getMass();
-    if(W < 2*mdec)
-    {
-       cout<<" ERROR: W="<<W<<endl;
-       iFbadevent = 1;
-       return;
-    }
-    pmag = sqrt(W*W/4. - mdec*mdec);
-
-    //     pick an orientation, based on the spin
-    //      phi has a flat distribution in 2*pi
-    phi = randyInstance.Rndom()*2.*starlightConstants::pi; //(random()/(RAND_MAX+1.0))* 2.*starlightConstants::pi;
-
-    //     find theta, the angle between one of the outgoing particles and
-    //    the beamline, in the frame of the two photons
-
-    if(getSpin()==0.5){
-       //  calculate a table of integrated angle values for leptons
-       //  JN05: Go from 1000->20000bins, use double precision for anglelep and thetalep. needed when W >6Gev.
-       hiresw = W;
-
-       anglelep[0] = 0.;
-
-       for(int i =1;i<=20000;i++)
-       {
-           hirestheta = starlightConstants::pi * double(i) /20000.;
-
-           //  Added sin(theta) phase space factor (not in thetalep) and changed E to W in thetalep call
-           //  11/9/2000 SRK
-           //  Note that thetalep is form invariant, so it could be called for E, theta_lab just
-           //  as well as W,theta_cm.  Since there is a boost from cm to lab below, the former is fine.
-
-           anglelep[i] = anglelep[i-1] + thetalep(hiresw,hirestheta)*sin(hirestheta);
-       }
-
-       hirestheta = 0.;
-       xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
-       hirestest = xtest;
-       for(int i =1;i<=20000;i++)
-       {
-           if(xtest > (anglelep[i]/anglelep[20000]))
-                hirestheta = starlightConstants::pi * double(i) / 20000.;
-       }
-       theta=hirestheta;
-
-    }
-
-    if(getSpin() != 0.5)
-        cout<<" This model cannot yet handle this spin value for lepton pairs: "<<getSpin()<<endl; 
-
-
-    //     compute unboosted momenta
-    px1 = sin(theta)*cos(phi)*pmag;
-    py1 = sin(theta)*sin(phi)*pmag;
-    pz1 = cos(theta)*pmag;
-    px2 = -px1;
-    py2 = -py1;
-    pz2 = -pz1;
-
-    //        compute energies
-    //Changed mass to W 11/9/2000 SRK
-    Ecm = sqrt(W*W+px0*px0+py0*py0+pz0*pz0);
-
-    E1 = sqrt(mdec*mdec+px1*px1+py1*py1+pz1*pz1);
-    E2 = sqrt(mdec*mdec+px2*px2+py2*py2+pz2*pz2);
-    //        decay tau to electrons
-    //        note that after this routine px1, etc., refer to the electrons
-    if(_GGlepInputpidtest == starlightConstants::TAUON)
-        tauDecay(px1,py1,pz1,E1,px2,py2,pz2,E2);
-
-    //     Lorentz transform into the lab frame
-    // betax,betay,betaz are the boost of the complete system
-    betax = -(px0/Ecm);
-    betay = -(py0/Ecm);
-    betaz = -(pz0/Ecm);
-
-    transform (betax,betay,betaz,E1,px1,py1,pz1,iFbadevent);
-    transform (betax,betay,betaz,E2,px2,py2,pz2,iFbadevent);
-
-
-    if(iFbadevent == 1)
-        return;
-
-    // change particle id from that of parent to that of daughters
-    // change taoun id into electron id, already switched particles in tau decay
-    if(_GGlepInputpidtest == starlightConstants::TAUON)
-        ipid = starlightConstants::ELECTRON;
-    //        electrons remain electrons; muons remain muons
-    if ((_GGlepInputpidtest == starlightConstants::ELECTRON) || (_GGlepInputpidtest == starlightConstants::MUON))
-        ipid = _GGlepInputpidtest;
-}
-
-
-//______________________________________________________________________________
-double Gammagammaleptonpair::thetalep(double W,double theta)
-{
-    //     This function calculates the cross-section as a function of
-    //     angle for a given W and Y, for the production of two muons.
-    //     (or tauons)
-    //    expression is taken from Brodsky et al. PRD 1971, 1532
-    //     equation 5.7
-    //     factor that are not dependant on theta are scrapped, so the
-    //     the absolute crosssections given by this function are inaccurate
-    //     here we are working in the CM frame of the photons and the last
-    //     term is 0
-
-    //    real function thetalep (W,theta)
-
-    double moverw=0., W1sq=0.;
-    double thetalep_r=0.,denominator=0.;
-
-    W1sq = (W / 2.)*(W/2.);
-    moverw = getMass()*getMass() / W1sq;
-    denominator = (1. - (1. - moverw)*(cos(theta)*cos(theta)));
-
-    thetalep_r = 2. + 4.*(1.-moverw)*((1.-moverw)*sin(theta)*sin(theta)*cos(theta)*cos(theta) + moverw) / (denominator*denominator);
-
-    return thetalep_r;
-
-}
-
-
-//______________________________________________________________________________
-starlightConstants::event Gammagammaleptonpair::produceEvent(int &ievent)
-{//returns the vector with the decay particles inside.
-    starlightConstants::event leptonpair; //This object will store all the tracks for a single event
-    double comenergy = 0.;
-    double rapidity = 0.;
-    double pairE = 0.;
-    double pairmomx=0.,pairmomy=0.,pairmomz=0.;
-    int iFbadevent=0;
-    starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
-       
-    double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.;
-//this function decays particles and writes events to a file
-    //zero out the event structure
-    leptonpair._numberOfTracks=0;
-    for(int i=0;i<4;i++)
-    {
-       leptonpair.px[i]=0.;
-       leptonpair.py[i]=0.;
-       leptonpair.pz[i]=0.;
-       leptonpair._fsParticle[i]=starlightConstants::UNKNOWN;
-       leptonpair._charge[i]=0;
-    }
-
-    pickw(comenergy);
-
-    picky(rapidity);
-
-    pairMomentum(comenergy,rapidity,pairE,pairmomx,pairmomy,pairmomz);
-    twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
-    if (iFbadevent==0){
-       int q1=0,q2=0; 
-
-       double xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
-       if (xtest<0.5)
-       {
-           q1=1;
-           q2=-1;
-       }
-       else{
-           q1=-1;
-           q2=1;
-       }       
-       leptonpair._numberOfTracks=2;//leptonpairs are two tracks...
-       leptonpair.px[0]=px1;
-       leptonpair.py[0]=py1;
-       leptonpair.pz[0]=pz1;
-       leptonpair._fsParticle[0]=ipid; 
-       leptonpair._charge[0]=q1;
-
-       leptonpair.px[1]=px2;
-       leptonpair.py[1]=py2;
-       leptonpair.pz[1]=pz2;
-       leptonpair._fsParticle[1]=ipid;
-       leptonpair._charge[1]=q2;
-
-       ievent=ievent+1;
-    }
-
-    return leptonpair;
-}
-
-
-//______________________________________________________________________________
-upcEvent Gammagammaleptonpair::produceEvent()
-{
-//returns the vector with the decay particles inside.
-
-   upcEvent event;
-
-   double comenergy = 0.;
-   double rapidity = 0.;
-   double pairE = 0.;
-   double pairmomx=0.,pairmomy=0.,pairmomz=0.;
-   int iFbadevent=0;
-   starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
-   
-   double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.,E2=0.,E1=0.;
-   bool accepted = false;
-   do{ 
-     //this function decays particles and writes events to a file
-     //zero out the event structure
-     pickw(comenergy);
-     
-     picky(rapidity);
-     
-     pairMomentum(comenergy,rapidity,pairE,pairmomx,pairmomy,pairmomz);
-   
-  
-     _nmbAttempts++;
-     twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
-     double pt1chk = sqrt(px1*px1+py1*py1);
-     double pt2chk = sqrt(px2*px2+py2*py2);
-    
-     double eta1 = pseudoRapidity(px1, py1, pz1);
-     double eta2 = pseudoRapidity(px2, py2, pz2);
-    
-     if(_ptCutEnabled && !_etaCutEnabled){
-       if(pt1chk > _ptCutMin && pt1chk < _ptCutMax &&  pt2chk > _ptCutMin && pt2chk < _ptCutMax){
-        accepted = true;
-        _nmbAccepted++;
-       }
-     }
-     else if(!_ptCutEnabled && _etaCutEnabled){
-       if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
-        accepted = true;
-        _nmbAccepted++;
-       }
-     }
-     else if(_ptCutEnabled && _etaCutEnabled){
-       if(pt1chk > _ptCutMin && pt1chk < _ptCutMax &&  pt2chk > _ptCutMin && pt2chk < _ptCutMax){
-        if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
-          accepted = true;
-           _nmbAccepted++;
-        }
-       }
-     }
-    
-   }while((_ptCutEnabled || _etaCutEnabled) && !accepted);
-   //twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
-   
-   if (iFbadevent==0){
-     int q1=0,q2=0; 
-     
-     double xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
-     if (xtest<0.5)
-       {
-        q1=1;
-        q2=-1;
-       }
-     else{
-       q1=-1;
-       q2=1;
-     }
-
-     // The new stuff
-     double mlepton = getMass(); 
-     E1 = sqrt( mlepton*mlepton + px1*px1 + py1*py1 + pz1*pz1 ); 
-     E2 = sqrt( mlepton*mlepton + px2*px2 + py2*py2 + pz2*pz2 ); 
-
-     starlightParticle particle1(px1, py1, pz1, E1, starlightConstants::UNKNOWN, -q1*ipid, q1);
-     event.addParticle(particle1);
-     
-     starlightParticle particle2(px2, py2, pz2, E2, starlightConstants::UNKNOWN, -q2*ipid, q2);
-     event.addParticle(particle2);
-     
-    }
-   return event;
-}
-
-
-//______________________________________________________________________________
-void Gammagammaleptonpair::calculateTable()
-{
-    //     this subroutine calculates the tables that are used
-    //     elsewhere in the montecarlo
-    //     the tauon decay is taken from V-A theory, 1 - 1/3 cos(theta)
-    //     the energy of each of the two leptons in tau decay
-    //     is calculated using formula 10.35 given
-    //     in introduction to elementary particles by D. griffiths
-    //     which assmes that the mass of the electron is 0.
-    //     the highest energy attainable by the electron in such a system is
-    //     .5 * mass of the tau
-
-    //    subroutine calculateTable
-
-
-    double E,theta;
-
-    _tautolangle[0] = 0.;
-    _dgammade[0] = 0.;
-
-
-    for(int i =1;i<=100;i++)
-    {
-       //     calculate energy of tau decay
-       E = double(i)/100. * .5 * starlightConstants::tauMass;
-       _dgammade[i] = _dgammade[i-1] + E*E * (1. - 4.*E/(3.*starlightConstants::tauMass));
-
-       //     calculate angles for tau
-       theta = starlightConstants::pi * double(i) / 100.;
-       _tautolangle[i] = _tautolangle[i-1] + (1 + 0.3333333 * cos(theta));
-    }
-
-
-}
-
-
-//______________________________________________________________________________
-void Gammagammaleptonpair::tauDecay(double &px1,double &py1,double &pz1,double &E1,double &px2,double &py2,double &pz2,double &E2)
-{
-    //     this routine assumes that the tauons decay to electrons and
-    //     calculates the directons of the decays
-
-    double Ee1,Ee2,theta1,theta2,phi1,phi2, ran1, ran2 ;
-    double pmag1,pex1,pey1,pez1,pex2,pey2,pez2,pmag2;
-    double betax,betay,betaz,dir;
-
-    int Tmp_Par=0; // temp variable for the transform function .. kind of starnge - being called with 7 parameter instead of 8
-
-    //     the highest energy attainable by the electron in this system is
-    //     .5 * mass of the tau
-
-    //     get two random numbers to compare with
-
-
-    ran1 = randyInstance.Rndom()*_dgammade[100];//(random()/(RAND_MAX+1.0)) * _dgammade[100];
-    ran2 = randyInstance.Rndom()*_dgammade[100];//(random()/(RAND_MAX+1.0)) * _dgammade[100];
-
-    //     compute the energies that correspond to those numbers
-    Ee1 = 0.;
-    Ee2 = 0.;
-
-    for( int i =1;i<=100;i++)
-    {
-       if (ran1 > _dgammade[i])
-           Ee1 = double(i) /100. * .5 * getMass();
-       if (ran2 > _dgammade[i])
-           Ee2 = double(i) /100. * .5 * getMass();
-    }
-
-    //     to find the direction of emmission, first
-    //     we determine if the tauons have spin of +1 or -1 along the
-    //     direction of the beam line
-    dir = 1.;
-    if ( randyInstance.Rndom() < 0.5 )//(random()/(RAND_MAX+1.0)) < 0.5)
-       dir = -1.;
-
-    //     get two random numbers to compare with
-    ran1 = randyInstance.Rndom()*_tautolangle[99];//(random()/(RAND_MAX+1.0))  * _tautolangle[100];
-    ran2 = randyInstance.Rndom()*_tautolangle[99];//(random()/(RAND_MAX+1.0))  * _tautolangle[100];
-
-    //     find the angles corrsponding to those numbers
-    theta1 = 0.;
-    theta2 = 0.;
-    for( int i =1;i<=100;i++)
-    {
-       if (ran1 > _tautolangle[i]) theta1 = starlightConstants::pi * double(i) /100.;
-       if (ran2 > _tautolangle[i]) theta2 = starlightConstants::pi * double(i) /100.;
-    }
-
-    //     grab another two random numbers to determine phi's
-    phi1 = randyInstance.Rndom()*2.*starlightConstants::pi;// (random()/(RAND_MAX+1.0))* 2. * starlightConstants::pi;
-    phi2 = randyInstance.Rndom()*2.*starlightConstants::pi;// (random()/(RAND_MAX+1.0))* 2. * starlightConstants::pi;
-    //     figure out the momenta of the electron in the frames of the
-    //     tauons from which they decayed, that is electron1 is in the
-    //     rest frame of tauon1 and e2 is in the rest fram of tau2
-    //     again the electrons are assumed to be massless
-    pmag1 = Ee1;
-    pex1 = cos(phi1)*sin(theta1)*pmag1;
-    pey1 = sin(phi1)*sin(theta1)*pmag1;
-    pez1 = cos(theta1)*pmag1*dir;
-    pmag2 = Ee2;
-    pex2 = cos(phi2)*sin(theta2)*pmag2;
-    pey2 = sin(phi2)*sin(theta2)*pmag2;
-    pez2 = cos(theta2)*pmag2*(-1.*dir);
-    //     now Lorentz transform into the frame of each of the particles
-    //     do particle one first
-    betax = -(px1/E1);
-    betay = -(py1/E1);
-    betaz = -(pz1/E1);
-//cout<<"2decay betax,pex1= "<<betax<<" "<<pex1<<endl;
-    transform (betax,betay,betaz,Ee1,pex1,pey1,pez1,Tmp_Par);
-    //     then do particle two
-    betax = -(px1/E2);
-    betay = -(py1/E2);
-    betaz = -(pz1/E2);
-
-    transform (betax,betay,betaz,Ee2,pex2,pey2,pez2, Tmp_Par);
-    //     finally dump the electron values into the approriate
-    //     variables
-    E1 = Ee1;
-    E2 = Ee2;
-    px1 = pex1;
-    px2 = pex2;
-    py1 = pey1;
-    py2 = pey2;
-    pz1 = pez1;
-    pz2 = pez2;
-}
-
-
-//______________________________________________________________________________
-double Gammagammaleptonpair::getMass()
-{
-    double leptonmass=0.;
-    switch(_GGlepInputpidtest){
-    case starlightConstants::ELECTRON:
-       leptonmass=starlightConstants::mel;
-       break;
-    case starlightConstants::MUON:
-       leptonmass=starlightConstants::muonMass;
-       break;
-    case starlightConstants::TAUON:
-       leptonmass=starlightConstants::tauMass;
-       break;
-    default:
-       cout<<"Not a recognized lepton, Gammagammaleptonpair::getmass(), mass = 0."<<endl;
-    }
-
-    return leptonmass;
-}
-
-
-//______________________________________________________________________________
-double Gammagammaleptonpair::getWidth()
-{
-    double leptonwidth=0.;
-    return leptonwidth;
-
-}
-
-
-//______________________________________________________________________________
-double Gammagammaleptonpair::getSpin()
-{
-    double leptonspin=0.5;
-    return leptonspin;
-}
diff --git a/STARLIGHT/starlight/src/inputParameters.cpp~ b/STARLIGHT/starlight/src/inputParameters.cpp~
deleted file mode 100644 (file)
index 607fbbd..0000000
+++ /dev/null
@@ -1,572 +0,0 @@
-///////////////////////////////////////////////////////////////////////////
-//
-//    Copyright 2010
-//
-//    This file is part of starlight.
-//
-//    starlight 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 3 of the License, or
-//    (at your option) any later version.
-//
-//    starlight 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 starlight. If not, see <http://www.gnu.org/licenses/>.
-//
-///////////////////////////////////////////////////////////////////////////
-//
-// File and Version Information:
-// $Rev:: 166                         $: revision of last commit
-// $Author:: odjuvsla                 $: author of last commit
-// $Date:: 2013-10-06 16:18:12 +0200 #$: date of last commit
-//
-// Description:
-//
-//
-//
-///////////////////////////////////////////////////////////////////////////
-
-
-#include <iostream>
-#include <fstream>
-
-#include "reportingUtils.h"
-#include "starlightconstants.h"
-#include "inputParameters.h"
-#include "inputParser.h"
-#include "starlightconfig.h"
-#include <cmath>
-#include <cstring>
-#include "randomgenerator.h"
-
-using namespace std;
-using namespace starlightConstants;
-
-parameterlist parameterbase::_parameters;
-
-#define REQUIRED true
-#define NOT_REQUIRED false
-
-//______________________________________________________________________________
-inputParameters::inputParameters()
-       : _configFileName        ("slight.in"),
-         _beam1Z                ("BEAM_1_Z",0),
-         _beam1A                ("BEAM_1_A",0),
-         _beam2Z                ("BEAM_2_Z",0),
-         _beam2A                ("BEAM_2_A",0),
-         _beam1LorentzGamma     ("BEAM_1_GAMMA",0),
-         _beam2LorentzGamma     ("BEAM_2_GAMMA",0),
-         _maxW                  ("W_MAX",0),
-         _minW                  ("W_MIN",0),
-         _nmbWBins              ("W_N_BINS",0),
-         _maxRapidity           ("RAP_MAX",0),
-         _nmbRapidityBins       ("RAP_N_BINS",0),
-         _ptCutEnabled          ("CUT_PT",false),
-         _ptCutMin              ("PT_MIN",0),
-         _ptCutMax              ("PT_MAX",0),
-         _etaCutEnabled         ("CUT_ETA",false),
-         _etaCutMin             ("ETA_MIN",0),
-         _etaCutMax             ("ETA_MAX",0),
-         _productionMode        ("PROD_MODE",0),
-         _nmbEventsTot          ("N_EVENTS",0),
-         _prodParticleId        ("PROD_PID",0),
-         _randomSeed            ("RND_SEED",0),
-         _outputFormat          ("OUTPUT_FORMAT",0),
-         _beamBreakupMode       ("BREAKUP_MODE",0),
-         _interferenceEnabled   ("INTERFERENCE",false),
-         _interferenceStrength  ("IF_STRENGTH",0),
-         _coherentProduction    ("COHERENT",false),
-         _incoherentFactor      ("INCO_FACTOR",0),
-         _deuteronSlopePar      ("BFORD",0),
-         _maxPtInterference     ("INT_PT_MAX",0),
-         _nmbPtBinsInterference ("INT_PT_N_BINS",0),
-         _ptBinWidthInterference("INT_PT_WIDTH",0),
-         _protonEnergy          ("PROTON_ENERGY",0, NOT_REQUIRED),
-         _minGammaEnergy        ("MIN_GAMMA_ENERGY",0, NOT_REQUIRED),
-         _maxGammaEnergy        ("MAX_GAMMA_ENERGY",0, NOT_REQUIRED),
-         _pythiaParams          ("PYTHIA_PARAMS","", NOT_REQUIRED),
-         _pythiaFullEventRecord ("PYTHIA_FULL_EVENTRECORD",false, NOT_REQUIRED),
-         _xsecCalcMethod        ("XSEC_METHOD",0, NOT_REQUIRED),
-         _nThreads              ("N_THREADS",1, NOT_REQUIRED),
-          _nBinsQKniehl          ("N_BINS_Q_KNIEHL", 0, NOT_REQUIRED),
-          _nBinsEKniehl          ("N_BINS_E_KNIEHL", 0, NOT_REQUIRED),
-          _nBinsBKniehl          ("N_BINS_B_KNIEHL", 0, NOT_REQUIRED),
-          _qMaxKniehl            ("Q_MAX_KNIEHL", 0, NOT_REQUIRED),
-          _eGammaMinKniehl       ("E_GAMMA_MIN_KNIEHL", 0, NOT_REQUIRED),
-          _eGammaMaxKniehl       ("E_GAMMA_MAX_KNIEHL", 0, NOT_REQUIRED),
-          _bMinKniehl            ("B_MIN_KNIEHL", 0, NOT_REQUIRED),
-          _bMaxKniehl            ("B_MAX_KNIEHL", 0, NOT_REQUIRED) 
-{
-  // All parameters must be initialised in initialisation list! 
-  // If not: error: 'parameter<T, validate>::parameter() [with T = unsigned int, bool validate = true]' is private
-  // or similar
-       
-       _ip.addParameter(_beam1Z);
-       _ip.addParameter(_beam2Z);
-       _ip.addParameter(_beam1A);
-       _ip.addParameter(_beam2A);
-
-       _ip.addParameter(_beam1LorentzGamma);
-       _ip.addParameter(_beam2LorentzGamma);
-       
-       _ip.addParameter(_maxW);
-       _ip.addParameter(_minW);
-       
-       _ip.addParameter(_nmbWBins);
-
-       _ip.addParameter(_maxRapidity);
-       _ip.addParameter(_nmbRapidityBins);
-       
-       _ip.addParameter(_ptCutEnabled);
-       _ip.addParameter(_ptCutMin);
-       _ip.addParameter(_ptCutMax);
-       
-       _ip.addParameter(_etaCutEnabled);
-       _ip.addParameter(_etaCutMax);
-       _ip.addParameter(_etaCutMin);
-       
-       _ip.addParameter(_productionMode);
-       _ip.addParameter(_nmbEventsTot);
-       _ip.addParameter(_prodParticleId);
-       _ip.addParameter(_randomSeed);
-       _ip.addParameter(_outputFormat);
-       _ip.addParameter(_beamBreakupMode);
-       _ip.addParameter(_interferenceEnabled);
-       _ip.addParameter(_interferenceStrength);
-       _ip.addParameter(_coherentProduction);
-       _ip.addParameter(_incoherentFactor);
-       _ip.addParameter(_deuteronSlopePar);
-       _ip.addParameter(_maxPtInterference);
-       _ip.addParameter(_nmbPtBinsInterference);
-       _ip.addParameter(_minGammaEnergy);
-       _ip.addParameter(_maxGammaEnergy);
-       _ip.addParameter(_pythiaParams);
-       _ip.addParameter(_pythiaFullEventRecord);
-       _ip.addParameter(_xsecCalcMethod);
-       _ip.addParameter(_nThreads);
-       _ip.addParameter(_nBinsBKniehl);
-       _ip.addParameter(_nBinsQKniehl);
-       _ip.addParameter(_nBinsEKniehl);
-       _ip.addParameter(_qMaxKniehl);
-       _ip.addParameter(_eGammaMaxKniehl);
-       _ip.addParameter(_eGammaMinKniehl);
-       _ip.addParameter(_bMaxKniehl);
-       _ip.addParameter(_bMinKniehl);
-}
-
-
-//______________________________________________________________________________
-inputParameters::~inputParameters()
-{ }
-
-
-//______________________________________________________________________________
-bool
-inputParameters::configureFromFile(const std::string &configFileName)
-{
-       // open config file
-       _configFileName = configFileName;
-       
-       int nParameters = _ip.parseFile(_configFileName);
-       
-       if(nParameters == -1) 
-       {
-               printWarn << "could not open file '" << _configFileName << "'" << endl;
-         return false;
-       }
-       
-       //ip.printParameterInfo(cout);
-       
-       if(_ip.validateParameters(cerr))
-               printInfo << "successfully read input parameters from '" << _configFileName << "'" << endl;
-       else {
-               printWarn << "problems reading input parameters from '" << _configFileName << "'" << endl
-                         << *this;
-               return false;
-       }
-       return true;
-}
- bool inputParameters::init()
- {
-       // Set the seed for the random generator
-       randyInstance.SetSeed(_randomSeed.value());
-       
-       // Calculate beam gamma in CMS frame
-       double rap1 = acosh(beam1LorentzGamma());
-       double rap2 = -acosh(beam2LorentzGamma());
-       _beamLorentzGamma = cosh((rap1-rap2)/2);
-       
-       std::cout << "Rapidity beam 1: " << rap1 << ", rapidity beam 2: " << rap2 << ", rapidity CMS system: " << (rap1+rap2)/2 << ", beam gamma in CMS: " << _beamLorentzGamma<< std::endl;
-       _ptBinWidthInterference = maxPtInterference() / nmbPtBinsInterference();
-       _protonEnergy           = _beamLorentzGamma * protonMass;
-
-       // define interaction type
-       switch (productionMode()) {
-       case 1:
-               _interactionType = PHOTONPHOTON;
-               break;
-       case 2:
-               _interactionType = PHOTONPOMERONNARROW;
-               break;
-       case 3:
-               _interactionType = PHOTONPOMERONWIDE;
-               break;
-        case 4:
-                _interactionType = PHOTONPOMERONINCOHERENT;
-                break;
-       case 5:
-               _interactionType = PHOTONUCLEARSINGLE;
-               break;
-       case 6:
-               _interactionType = PHOTONUCLEARDOUBLE;
-               break;
-       case 7:
-               _interactionType = PHOTONUCLEARSINGLEPA;
-               break;
-       case 8:
-               _interactionType = PHOTONUCLEARSINGLEPAPY;
-               break;
-//     case 9:
-//             _interactionType = PHOTONPHOTONKNIEHL;
-//             break;
-//         case 10:
-//                 _interactionType = PHOTONPHOTONKNIEHLMODIFIED;
-//                 break;
-       default:
-               printWarn << "unknown production mode '" << _productionMode << "'" << endl;
-               return false;
-       }
-  
-       //Trying to define the proper Wmins and Wmaxs. a TEMPORARY fix....Better solution=??
-       double mass        = 0;
-       double width       = 0;
-       double defaultMinW = 0;  // default for _minW, unless it is defined later [GeV/c^2]
-       switch (prodParticleId()) {
-       case 11:  // e+e- pair
-               _particleType = ELECTRON;
-               _decayType    = LEPTONPAIR;
-               defaultMinW   = 0.01; // default is 0.01; up to 0.15 is safe for Summer 2000 triggering for e+e- pairs
-               break;
-       case 13:  // mu+mu- pair
-               _particleType = MUON;
-               _decayType    = LEPTONPAIR;
-               defaultMinW   = 2 * muonMass;
-               break;
-       case 15:  // tau+tau- pair
-               _particleType = TAUON;
-               _decayType    = LEPTONPAIR;
-               defaultMinW   = 2 * tauMass;
-               break;
-//     case 24:  // W+W- pair
-//             _particleType = W;
-//             _decayType    = WW;
-//             defaultMinW   = 2 * muonMass;
-//             break;
-       case 115:  // a_2(1320)
-               _particleType = A2;
-               _decayType    = SINGLEMESON;
-               break;
-       case 221:  // eta
-               _particleType = ETA;
-               _decayType    = SINGLEMESON;
-               break;
-       case 225:  // f_2(1270)
-               _particleType = F2;
-                defaultMinW   = 2*pionChargedMass;
-               _decayType    = SINGLEMESON;
-               break;
-       case 331:  // eta'(958)
-               _particleType = ETAPRIME;
-               _decayType    = SINGLEMESON;
-               break;
-       case 335:  // f_2'(1525)
-               _particleType = F2PRIME;
-               _decayType    = SINGLEMESON;
-               break;
-       case 441:  // eta_c(1s)
-               _particleType = ETAC;
-               _decayType    = SINGLEMESON;
-                defaultMinW   = etaCMass - 5 * 0.0267;
-               break;
-       case 9010221:  // f_0(980), was orginally called 10221? updated to standard number
-               _particleType = F0;
-               _decayType    = SINGLEMESON;
-                defaultMinW   = 2*pionNeutralMass;
-               break;
-       case 33:  // Z"/Z03
-               _particleType = ZOVERZ03;
-               _decayType    = SINGLEMESON;
-               break;
-//     case 25: // Higgs
-//             _particleType = HIGGS;
-//             _decayType    = SINGLEMESON;
-//             break;
-       case 113:  // rho(770)
-               _particleType = RHO;
-               _decayType    = WIDEVMDEFAULT;
-               mass          = 0.7685;
-               width         = 0.1507;
-               defaultMinW   = 2 * pionChargedMass;
-               _maxW         = mass + 5 * width;
-               break;
-       case 913:  // rho(770) with direct pi+pi- decay, interference given by ZEUS data
-               _particleType = RHOZEUS;
-               _decayType    = WIDEVMDEFAULT;
-               mass          = 0.7685;                    
-               width         = 0.1507;
-               defaultMinW   = 2 * pionChargedMass;
-               _maxW         = mass + 5 * width;  // use the same 1.5GeV max mass as ZEUS
-               break;
-       case 999:  // pi+pi-pi+pi- phase space decay
-               _particleType = FOURPRONG;
-               _decayType    = WIDEVMDEFAULT;
-               mass          = 1.350;
-               width         = 0.360;
-               defaultMinW   = 4 * pionChargedMass;
-               _maxW         = 3;
-               break;
-       case 223:  // omega(782)
-               _particleType = OMEGA;
-               _decayType    = NARROWVMDEFAULT;  // will probably be moved to 3-body decay
-               mass          = 0.78194;
-               width         = 0.00843;
-               defaultMinW   = mass - 5 * width;
-               _maxW         = mass + 5 * width;
-               break;
-       case 333:  // phi(1020)
-               _particleType = PHI;
-               _decayType    = NARROWVMDEFAULT;
-               mass          = 1.019413;
-               width         = 0.00443;
-               defaultMinW   = 2 * kaonChargedMass;
-               _maxW         = mass + 5 * width;
-               break;
-       case 443:  // J/psi
-               _particleType = JPSI;
-               _decayType    = NARROWVMDEFAULT;
-               mass          = 3.09692;   // JN  3.09688;
-               width         = 0.000091;  // JN  0.000087;
-               defaultMinW   = mass - 5 * width;
-               _maxW         = mass + 5 * width;
-               break;
-       case 443011:  // J/psi
-               _particleType = JPSI_ee;
-               _decayType    = NARROWVMDEFAULT;
-               mass          = 3.09692;   // JN  3.09688;
-               width         = 0.000091;  // JN  0.000087;
-               defaultMinW   = mass - 5 * width;
-               _maxW         = mass + 5 * width;
-               break;
-       case 443013:  // J/psi
-               _particleType = JPSI_mumu;
-               _decayType    = NARROWVMDEFAULT;
-               mass          = 3.09692;   // JN  3.09688;
-               width         = 0.000091;  // JN  0.000087;
-               defaultMinW   = mass - 5 * width;
-               _maxW         = mass + 5 * width;
-               break;
-       case 444:  // J/psi
-               _particleType = JPSI2S;
-               _decayType    = NARROWVMDEFAULT;
-               mass          = 3.686093;
-               width         = 0.000337;
-               defaultMinW   = mass - 5 * width;
-               _maxW         = mass + 5 * width;
-               break;
-       case 444011: // J/psi
-               _particleType = JPSI2S_ee;
-               _decayType    = NARROWVMDEFAULT;
-               mass          = 3.686093;     
-               width         = 0.000337;
-               defaultMinW   = mass - 5 * width;
-               _maxW         = mass + 5 * width;
-               break;
-       case 444013:  // J/psi
-               _particleType = JPSI2S_mumu;
-               _decayType    = NARROWVMDEFAULT;
-               mass          = 3.686093;
-               width         = 0.000337;
-               defaultMinW   = mass - 5 * width;
-               _maxW         = mass + 5 * width;
-               break;
-       case 553:  // Upsilon
-               _particleType = UPSILON;
-               _decayType    = NARROWVMDEFAULT;
-               mass          = 9.46030;
-               width         = 0.00005402;
-               defaultMinW   = mass - 5 * width;
-               _maxW         = mass + 5 * width;
-               break;
-       case 553011:  // Upsilon
-               _particleType = UPSILON_ee;
-               _decayType    = NARROWVMDEFAULT;
-               mass          = 9.46030;
-               width         = 0.00005402;
-               defaultMinW   = mass - 5 * width;
-               _maxW         = mass + 5 * width;
-               break;
-       case 553013:  // Upsilon
-               _particleType = UPSILON_mumu;
-               _decayType    = NARROWVMDEFAULT;
-               mass          = 9.46030;
-               width         = 0.00005402;
-               defaultMinW   = mass - 5 * width;
-               _maxW         = mass + 5 * width;
-               break;
-       case 554:  // Upsilon(2S)
-               _particleType = UPSILON2S;
-               _decayType    = NARROWVMDEFAULT;
-               mass          = 10.02326;
-               width         = 0.00003198;
-               defaultMinW   = mass - 5 * width;
-               _maxW         = mass + 5 * width;
-               break;
-       case 554011:  // Upsilon(2S)
-               _particleType = UPSILON2S_ee;
-               _decayType    = NARROWVMDEFAULT;
-               mass          = 10.02326;
-               width         = 0.00003198;
-               defaultMinW   = mass - 5 * width;
-               _maxW         = mass + 5 * width;
-               break;
-       case 554013:  // Upsilon(2S)
-               _particleType = UPSILON2S_mumu;
-               _decayType    = NARROWVMDEFAULT;
-               mass          = 10.02326;
-               width         = 0.00003198;
-               defaultMinW   = mass - 5 * width;
-               _maxW         = mass + 5 * width;
-               break;
-       case 555:  // Upsilon(3S)
-               mass          = 10.3552;
-               width         = 0.00002032;
-               defaultMinW   = mass - 5 * width;
-               _maxW         = mass + 5 * width;
-               _particleType = UPSILON3S;
-               _decayType    = NARROWVMDEFAULT;
-               break;
-       case 555011:  // Upsilon(3S)
-               _particleType = UPSILON3S_ee;
-               _decayType    = NARROWVMDEFAULT;
-               mass          = 10.3552;
-               width         = 0.00002032;
-               defaultMinW   = mass - 5 * width;
-               _maxW         = mass + 5 * width;
-               break;
-       case 555013:  // Upsilon(3S)
-               _particleType = UPSILON3S_mumu;
-               _decayType    = NARROWVMDEFAULT;
-               mass          = 10.3552;
-               width         = 0.00002032;
-               defaultMinW   = mass - 5 * width;
-               _maxW         = mass + 5 * width;
-               break;
-       default:
-               printWarn << "unknown particle ID " << _prodParticleId << endl;
-               return false;
-       }  // _prodParticleId
-
-       if (_minW.value() == -1)
-               _minW = defaultMinW;
-
-       printInfo << "using the following " << *this;
-       
-       return true;
-}
-
-
-//______________________________________________________________________________
-ostream&
-inputParameters::print(ostream& out) const
-{
-       out << "starlight parameters:" << endl
-           << "    config file name  ...................... '" << _configFileName << "'" << endl
-           << "    beam 1 atomic number ................... " << _beam1Z.value() << endl
-           << "    beam 1 atomic mass number .............. " << _beam1A.value() << endl
-           << "    beam 2 atomic number ................... " << _beam2Z.value() << endl
-           << "    beam 2 atomic mass number .............. " << _beam2A.value() << endl
-           << "    Lorentz gamma of beams in CM frame ..... " << _beamLorentzGamma << endl
-           << "    mass W of produced hadronic system ..... " << _minW.value() << " < W < " << _maxW.value() << " GeV/c^2" << endl
-           << "    # of W bins ............................ " << _nmbWBins.value() << endl
-           << "    maximum absolute value for rapidity .... " << _maxRapidity.value() << endl
-           << "    # of rapidity bins ..................... " << _nmbRapidityBins.value() << endl
-           << "    cut in pT............................... " << yesNo(_ptCutEnabled.value()) << endl
-           << "        minumum pT.......................... " << _ptCutMin.value() << " GeV/c" << endl
-           << "        maximum pT.......................... " << _ptCutMax.value() << " GeV/c" << endl
-           << "    cut in eta.............................. " << yesNo(_etaCutEnabled.value()) << endl
-           << "        minumum eta......................... " << _etaCutMin.value() << endl
-           << "        maximum eta......................... " << _etaCutMax.value() << endl
-           << "    meson production mode .................. " << _productionMode.value() << endl
-           << "    number of events to generate ........... " << _nmbEventsTot.value() << endl
-           << "    PDG ID of produced particle ............ " << _prodParticleId.value() << endl
-           << "    seed for random generator .............. " << _randomSeed.value() << endl
-           << "    output format .......................... " << _outputFormat.value() << endl
-           << "    breakup mode for beam particles ........ " << _beamBreakupMode.value() << endl
-           << "    interference enabled ................... " << yesNo(_interferenceEnabled.value()) << endl
-           << "    interference strength .................. " << _interferenceStrength.value() << endl
-           << "    coherent scattering off nucleus ........ " << yesNo(_coherentProduction.value()) << endl
-           << "    scaling factor for incoh. VM prod. ..... " << _incoherentFactor.value() << endl
-           << "    deuteron slope parameter ............... " << _deuteronSlopePar.value() << " (GeV/c)^{-2}" << endl
-           << "    maximum p_T for interference calc. ..... " << _maxPtInterference.value() << " GeV/c" << endl
-           << "    # of p_T bins for interference calc. ... " << _nmbPtBinsInterference.value() << endl;
-       return out;
-}
-
-
-//______________________________________________________________________________
-ostream&
-inputParameters::write(ostream& out) const
-{
-  
-       out << "BEAM_1_Z"      << beam1Z               () <<endl
-           << "BEAM_2_Z"      << beam1A               () <<endl
-           << "BEAM_1_A"      << beam2Z               () <<endl
-           << "BEAM_2_A"      << beam2A               () <<endl
-           << "BEAM_GAMMA"    << beamLorentzGamma     () <<endl
-           << "W_MAX"         << maxW                 () <<endl
-           << "W_MIN"         << minW                 () <<endl
-           << "W_N_BINS"      << nmbWBins             () <<endl
-           << "RAP_MAX"       << maxRapidity          () <<endl
-           << "RAP_N_BINS"    << nmbRapidityBins      () <<endl
-           << "CUT_PT"        << ptCutEnabled         () <<endl
-           << "PT_MIN"        << ptCutMin             () <<endl
-           << "PT_MAX"        << ptCutMax             () <<endl
-           << "CUT_ETA"       << etaCutEnabled        () <<endl
-           << "ETA_MIN"       << etaCutMin            () <<endl
-           << "ETA_MAX"       << etaCutMax            () <<endl
-           << "PROD_MODE"     << productionMode       () <<endl
-           << "N_EVENTS"      << nmbEvents            () <<endl
-           << "PROD_PID"      << prodParticleId       () <<endl
-           << "RND_SEED"      << randomSeed           () <<endl
-           << "OUTPUT_FORMAT" << outputFormat         () <<endl
-           << "BREAKUP_MODE"  << beamBreakupMode      () <<endl
-           << "INTERFERENCE"  << interferenceEnabled  () <<endl
-           << "IF_STRENGTH"   << interferenceStrength () <<endl
-           << "COHERENT"      << coherentProduction   () <<endl
-           << "INCO_FACTOR"   << incoherentFactor     () <<endl
-           << "BFORD"         << deuteronSlopePar     () <<endl
-           << "INT_PT_MAX"    << maxPtInterference    () <<endl
-           << "INT_PT_N_BINS" << nmbPtBinsInterference() <<endl;
-
-       return out;
-}
-
-std::string
-inputParameters::parameterValueKey() const 
-{
- // std::stringstream s;
-  
-//   s <<_beam1A<<_beam1Z<<_beam2A<<_beam1LorentzGamma<<_beam2LorentzGamma<<_maxW<<_minW;
-//     <<_nmbWBins<<_maxRapidity<<_nmbRapidityBins<<_
-    
-  
-  //return s;
-  
-  return parameterbase::_parameters.validationKey()
-  ;
-}
diff --git a/STARLIGHT/starlight/src/inputParser.cpp~ b/STARLIGHT/starlight/src/inputParser.cpp~
deleted file mode 100644 (file)
index 698fc6a..0000000
+++ /dev/null
@@ -1,341 +0,0 @@
-///////////////////////////////////////////////////////////////////////////
-//
-//    Copyright 2011
-//
-//    This file is part of starlight.
-//
-//    starlight 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 3 of the License, or
-//    (at your option) any later version.
-//
-//    starlight 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 starlight. If not, see <http://www.gnu.org/licenses/>.
-//
-///////////////////////////////////////////////////////////////////////////
-//
-// File and Version Information:
-// $Rev:: 155                       $: revision of last commit
-// $Author:: odjuvsla           $: author of last commit
-// $Date:: 2013-10-06 16:17:50 +#$: date of last commit
-//
-// Description:
-//
-//
-//
-///////////////////////////////////////////////////////////////////////////
-
-
-#include "../include/inputParser.h"
-#include <fstream>
-#include <cstdlib>
-#include <algorithm>
-
-inputParser::inputParser()
-{
-}
-
-inputParser::~inputParser()
-{
-}
-
-int inputParser::parseFile(std::string filename)
-{
-
-    std::ifstream infile(filename.c_str());
-    if ((!infile) || (!infile.good())) 
-    {
-      return -1;
-    }
-    
-    int lineSize = 256;
-    char tmp[lineSize];
-    int nParameters = 0;
-    while (!infile.getline(tmp, lineSize).eof())
-    {
-
-        std::string line(tmp);
-       nParameters += parseString(line);
-    }
-
-    infile.close();
-    return nParameters;
-
-}
-int inputParser::parseString(std::string str)
-{
-
-    std::string word;
-    std::string name;
-    std::string val;
-
-    std::map<std::string, _parameter<int> >::iterator intIt;
-    std::map<std::string, _parameter<unsigned int> >::iterator uIntIt;
-    std::map<std::string, _parameter<float> >::iterator floatIt;
-    std::map<std::string, _parameter<double> >::iterator doubleIt;
-    std::map<std::string, _parameter<bool> >::iterator boolIt;
-    std::map<std::string, _parameter<std::string> >::iterator stringIt;
-    
-    // Check if there is commented out stuff...
-    size_t pos = str.find_first_of("#");
-
-    // Cut the comment out of the str
-    if (pos != str.npos) str.erase(pos, str.find_first_of('\n'));
-
-    // Find the required equal sign and split the string into name and value
-    size_t eqPos = str.find("=");
-    std::string whitespaces (" \t\f\v\n\r");
-
-    name = "";
-    val = "";
-
-    if (eqPos != str.npos)
-    {
-        name = str.substr(0, eqPos);
-        name.erase(name.find_last_not_of(whitespaces)+1);
-        val = str.substr(eqPos+1);
-        val.erase(0, val.find_first_not_of(whitespaces));
-    }
-
-    if (name.length() > 0 && val.length() > 0)
-    {
-        intIt = _intParameters.find(name);
-        if (intIt != _intParameters.end())
-        {
-            intIt->second._found = true;
-            *(intIt->second._val) = atoi(val.c_str());
-           return true;
-        }
-        uIntIt = _uintParameters.find(name);
-        if (uIntIt != _uintParameters.end())
-        {
-            uIntIt->second._found = true;
-            *(uIntIt->second._val) = atoi(val.c_str());
-           return true;
-        }
-        floatIt = _floatParameters.find(name);
-        if (floatIt != _floatParameters.end())
-        {
-            floatIt->second._found = true;
-            *(floatIt->second._val) = atof(val.c_str());
-           return true;
-        }
-        doubleIt = _doubleParameters.find(name);
-        if (doubleIt != _doubleParameters.end())
-        {
-            doubleIt->second._found = true;
-            *(doubleIt->second._val) = atof(val.c_str());
-           return true;
-        }
-        boolIt = _boolParameters.find(name);
-        if (boolIt != _boolParameters.end())
-        {
-            boolIt->second._found = true;
-            *(boolIt->second._val) = atoi(val.c_str());
-           return true;
-        }
-        stringIt = _stringParameters.find(name);
-        if (stringIt != _stringParameters.end())
-        {
-            stringIt->second._found = true;
-            *(stringIt->second._val) = val;
-           return true;
-        }
-
-    }
-    return false;
-}
-void inputParser::addIntParameter(std::string name, int *var, bool required)
-{
-    _parameter<int> par(name, var, required);
-    _intParameters.insert(std::pair<std::string, _parameter<int> >(name, par));
-}
-
-void inputParser::addUintParameter(std::string name, unsigned int *var, bool required)
-{
-    _parameter<unsigned int> par(name, var, required);
-    _uintParameters.insert(std::pair<std::string, _parameter<unsigned int> >(name, par));
-}
-
-void inputParser::addFloatParameter(std::string name, float *var, bool required)
-{
-    _parameter<float> par(name, var, required);
-    _floatParameters.insert(std::pair<std::string, _parameter<float> >(name, par));
-}
-
-void inputParser::addDoubleParameter(std::string name, double *var, bool required)
-{
-    _parameter<double> par(name, var, required);
-    _doubleParameters.insert(std::pair<std::string, _parameter<double> >(name, par));
-}
-
-void inputParser::addBoolParameter(std::string name, bool *var, bool required)
-{
-    _parameter<bool> par(name, var, required);
-    _boolParameters.insert(std::pair<std::string, _parameter<bool> >(name, par));
-}
-
-void inputParser::addStringParameter(std::string name, std::string *var, bool required)
-{
-    _parameter<std::string> par(name, var, required);
-    _stringParameters.insert(std::pair<std::string, _parameter<std::string> >(name, par));
-}
-
-void inputParser::printParameterInfo(std::ostream &out)
-{
-
-    std::map<std::string, _parameter<int> >::iterator intIt;
-    std::map<std::string, _parameter<unsigned int> >::iterator uIntIt;
-    std::map<std::string, _parameter<float> >::iterator floatIt;
-    std::map<std::string, _parameter<double> >::iterator doubleIt;
-    std::map<std::string, _parameter<bool> >::iterator boolIt;
-    out << "#########################################" << std::endl;
-    out << "PARAMETER:\t\tVALUE:" << std::endl;
-    out << "#########################################" << std::endl;
-    out << "-----------------------------------------" << std::endl;
-    for (intIt = _intParameters.begin(); intIt != _intParameters.end(); ++intIt)
-    {
-        intIt->second.printParameterInfo();
-        out << "-----------------------------------------" << std::endl;
-    }
-    for (uIntIt = _uintParameters.begin(); uIntIt != _uintParameters.end(); ++uIntIt)
-    {
-        uIntIt->second.printParameterInfo();
-        out << "-----------------------------------------" << std::endl;
-    }
-    for (floatIt = _floatParameters.begin(); floatIt != _floatParameters.end(); ++floatIt)
-    {
-        floatIt->second.printParameterInfo();
-        out << "-----------------------------------------" << std::endl;
-    }
-    for (doubleIt = _doubleParameters.begin(); doubleIt != _doubleParameters.end(); ++doubleIt)
-    {
-        doubleIt->second.printParameterInfo();
-        out << "-----------------------------------------" << std::endl;
-    }
-    for (boolIt = _boolParameters.begin(); boolIt != _boolParameters.end(); ++boolIt)
-    {
-        boolIt->second.printParameterInfo();
-        out << "-----------------------------------------" << std::endl;
-    }
-    out << "#########################################" << std::endl;
-}
-
-bool inputParser::validateParameters(std::ostream& warnOut, std::ostream& errOut)
-{
-
-    int nNonCriticalMissing = 0;
-    int nCriticalMissing = 0;
-
-    std::map<std::string, _parameter<int> >::iterator intIt;
-    std::map<std::string, _parameter<float> >::iterator floatIt;
-    std::map<std::string, _parameter<double> >::iterator doubleIt;
-    std::map<std::string, _parameter<bool> >::iterator boolIt;
-
-    for (intIt = _intParameters.begin(); intIt != _intParameters.end(); ++intIt)
-    {
-        if (!intIt->second._found)
-        {
-            if (intIt->second._required)
-            {
-                errOut << "Could not find parameter: " << intIt->second._name << " which is required. Please specify this parameter in the config file!" << std::endl;
-                nCriticalMissing++;
-            }
-            else
-            {
-                warnOut << "Could not find parameter: " << intIt->second._name << ", but it is not required, using default value: " << *intIt->second._val << std::endl;
-                nNonCriticalMissing++;
-            }
-        }
-    }
-    for (floatIt = _floatParameters.begin(); floatIt != _floatParameters.end(); ++floatIt)
-    {
-        if (!floatIt->second._found)
-        {
-            if (floatIt->second._required)
-            {
-                errOut << "Could not find parameter: " << floatIt->second._name << " which is required. Please specify this parameter in the config file!" << std::endl;
-                nCriticalMissing++;
-            }
-            else
-            {
-                warnOut << "Could not find parameter: " << floatIt->second._name << ", but it is not required, using default value: " << *floatIt->second._val << std::endl;
-                nNonCriticalMissing++;
-            }
-        }
-    }
-    for (doubleIt = _doubleParameters.begin(); doubleIt != _doubleParameters.end(); ++doubleIt)
-    {
-        if (!doubleIt->second._found)
-        {
-            if (doubleIt->second._required)
-            {
-                errOut << "Could not find parameter: " << doubleIt->second._name << " which is required. Please specify this parameter in the config file!" << std::endl;
-                nCriticalMissing++;
-            }
-            else
-            {
-                warnOut << "Could not find parameter: " << doubleIt->second._name << ", but it is not required, using default value: " << *doubleIt->second._val << std::endl;
-                nNonCriticalMissing++;
-            }
-        }
-    }
-    for (boolIt = _boolParameters.begin(); boolIt != _boolParameters.end(); ++boolIt)
-    {
-        if (!boolIt->second._found)
-        {
-            if (boolIt->second._required)
-            {
-                errOut << "Could not find parameter: " << boolIt->second._name << " which is required. Please specify this parameter in the config file!" << std::endl;
-                nCriticalMissing++;
-            }
-            else
-            {
-                warnOut << "Could not find parameter: " << boolIt->second._name << ", but it is not required, using default value: " << *boolIt->second._val << std::endl;
-                nNonCriticalMissing++;
-            }
-        }
-    }
-    if(nCriticalMissing > 0) return false;
-    return true;
-}
-
-
-template<>
-void inputParser::addParameter(const std::string& name, int * varPtr, bool required)
-{
-  addIntParameter(name, varPtr, required);
-}
-template<>
-void inputParser::addParameter(const std::string& name, unsigned int * varPtr, bool required)
-{
-  addUintParameter(name, varPtr, required);
-}
-template<>
-void inputParser::addParameter(const std::string& name, float * varPtr, bool required)
-{
-  addFloatParameter(name, varPtr, required);
-}
-
-template<>
-void inputParser::addParameter(const std::string& name, double * varPtr, bool required)
-{
-  addDoubleParameter(name, varPtr, required);
-}
-
-template<>
-void inputParser::addParameter(const std::string& name, bool * varPtr, bool required)
-{
-  addBoolParameter(name, varPtr, required);
-}
-
-template<>
-void inputParser::addParameter(const std::string& name, std::string * varPtr, bool required)
-{
-  addStringParameter(name, varPtr, required);
-}
\ No newline at end of file
diff --git a/STARLIGHT/starlight/src/nBodyPhaseSpaceGen.cpp~ b/STARLIGHT/starlight/src/nBodyPhaseSpaceGen.cpp~
deleted file mode 100644 (file)
index 228d2d1..0000000
+++ /dev/null
@@ -1,258 +0,0 @@
-///////////////////////////////////////////////////////////////////////////
-//
-//    Copyright 2010
-//
-//    This file is part of starlight.
-//
-//    starlight 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 3 of the License, or
-//    (at your option) any later version.
-//       
-//    starlight 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 starlight. If not, see <http://www.gnu.org/licenses/>.
-//
-///////////////////////////////////////////////////////////////////////////
-//
-// File and Version Information:
-// $Rev:: 159                         $: revision of last commit
-// $Author:: odjuvsla                 $: author of last commit
-// $Date:: 2013-10-06 16:17:58 +0200 #$: date of last commit
-//
-// Description:
-//     see nBodyPhaseSpaceGen.h
-//
-//
-///////////////////////////////////////////////////////////////////////////
-
-
-#include <algorithm>
-
-#include "nBodyPhaseSpaceGen.h"
-
-
-using namespace std;
-using namespace starlightConstants;
-
-
-nBodyPhaseSpaceGen::nBodyPhaseSpaceGen()
-       : _n                (0),
-         _norm             (0),
-         _weight           (0),
-         _maxWeightObserved(0),
-         _maxWeight        (0)
-{ }
-
-
-nBodyPhaseSpaceGen::~nBodyPhaseSpaceGen()
-{ }
-
-
-// sets decay constants and prepares internal variables
-bool
-nBodyPhaseSpaceGen::setDecay(const vector<double>& daughterMasses)  // array of daughter particle masses
-{
-       _n = daughterMasses.size();
-       if (_n < 2) {
-               printWarn << "number of daughters = " << _n << " does not make sense." << endl;
-               return false;
-       }
-       // copy daughter masses
-       _m.clear();
-       _m = daughterMasses;
-       // prepare effective mass vector
-       _M.clear();
-       _M.resize(_n, 0);
-       _M[0] = _m[0];
-       // prepare angle vectors
-       _cosTheta.clear();
-       _cosTheta.resize(_n, 0);
-       _phi.clear();
-       _phi.resize(_n, 0);
-       // calculate daughter mass sums
-       _mSum.clear();
-       _mSum.resize(_n, 0);
-       _mSum[0] = _m[0];
-       for (unsigned int i = 1; i < _n; ++i)
-               _mSum[i] = _mSum[i - 1] + _m[i];
-       // prepare breakup momentum vector
-       _breakupMom.clear();
-       _breakupMom.resize(_n, 0);
-       // prepare vector for daughter Lorentz vectors
-       _daughters.clear();
-       _daughters.resize(_n, lorentzVector(0, 0, 0, 0));
-       // calculate normalization
-       _norm = 1 / (2 * pow(twoPi, 2 * (int)_n - 3) * factorial(_n - 2));
-       resetMaxWeightObserved();
-       return true;
-}
-
-
-// set decay constants and prepare internal variables
-bool
-nBodyPhaseSpaceGen::setDecay(const unsigned int nmbOfDaughters,  // number of daughter particles
-                             const double*      daughterMasses)  // array of daughter particle masses
-{
-       vector <double> m;
-       m.resize(nmbOfDaughters, 0);
-       for (unsigned int i = 0; i < nmbOfDaughters; ++i)
-               m[i] = daughterMasses[i];
-       return setDecay(m);
-}
-
-
-// generates event with certain n-body mass and momentum and returns event weigth
-// general purpose function
-double
-nBodyPhaseSpaceGen::generateDecay(const lorentzVector& nBody)  // Lorentz vector of n-body system in lab frame
-{
-       const double nBodyMass = nBody.M();
-       if (_n < 2) {
-               printWarn << "number of daughter particles = " << _n << " is smaller than 2. "
-                         << "weight is set to 0." << endl;
-               _weight = 0;
-       } else if (nBodyMass < _mSum[_n - 1]) {
-               printWarn << "n-body mass = " << nBodyMass << " is smaller than sum of daughter masses = "
-                         << _mSum[_n - 1] << ". weight is set to 0." << endl;
-               _weight = 0;
-       } else {
-               pickMasses(nBodyMass);
-               calcWeight();
-               pickAngles();
-               calcEventKinematics(nBody);
-       }
-       return _weight;
-}
-
-
-// generates full event with certain n-body mass and momentum only, when event is accepted (return value = true)
-// this function is more efficient, if only weighted evens are needed
-bool
-nBodyPhaseSpaceGen::generateDecayAccepted(const lorentzVector& nBody,      // Lorentz vector of n-body system in lab frame
-                                          const double         maxWeight)  // if positive, given value is used as maximum weight, otherwise _maxWeight
-{
-       const double nBodyMass = nBody.M();
-       if (_n < 2) {
-               printWarn << "number of daughter particles = " << _n << " is smaller than 2. "
-                         << "no event generated." << endl;
-               return false;
-       } else if (nBodyMass < _mSum[_n - 1]) {
-               printWarn << "n-body mass = " << nBodyMass << " is smaller than sum of daughter masses = "
-                         << _mSum[_n - 1] << ". no event generated." << endl;
-               return false;
-       }
-       pickMasses(nBodyMass);
-       calcWeight();
-       if (!eventAccepted(maxWeight))
-               return false;
-       pickAngles();
-       calcEventKinematics(nBody);
-       return true;
-}
-
-
-// randomly choses the (n - 2) effective masses of the respective (i + 1)-body systems
-void
-nBodyPhaseSpaceGen::pickMasses(const double nBodyMass)  // total energy of the system in its RF
-{
-       _M[_n - 1] = nBodyMass;
-       // create vector of sorted random values
-       vector<double> r(_n - 2, 0);  // (n - 2) values needed for 2- through (n - 1)-body systems
-       for (unsigned int i = 0; i < (_n - 2); ++i)
-               r[i] = random();
-       sort(r.begin(), r.end());
-       // set effective masses of (intermediate) two-body decays
-       const double massInterval = nBodyMass - _mSum[_n - 1];  // kinematically allowed mass interval
-       for (unsigned int i = 1; i < (_n - 1); ++i)             // loop over intermediate 2- to (n - 1)-bodies
-               _M[i] = _mSum[i] + r[i - 1] * massInterval;           // _mSum[i] is minimum effective mass
-}
-
-
-// computes event weight (= integrand value) and breakup momenta
-// uses vector of intermediate two-body masses prepared by pickMasses()
-double
-nBodyPhaseSpaceGen::calcWeight()
-{
-       for (unsigned int i = 1; i < _n; ++i)  // loop over 2- to n-bodies
-               _breakupMom[i] = breakupMomentum(_M[i], _M[i - 1], _m[i]);
-       double momProd = 1;                    // product of breakup momenta
-       for (unsigned int i = 1; i < _n; ++i)  // loop over 2- to n-bodies
-               momProd *= _breakupMom[i];
-       const double massInterval = _M[_n - 1] - _mSum[_n - 1];  // kinematically allowed mass interval
-       _weight = _norm * pow(massInterval, (int)_n - 2) * momProd / _M[_n - 1];
-       if (_weight > _maxWeightObserved)
-               _maxWeightObserved = _weight;
-       if (std::isnan(_weight))
-               printWarn << "weight = " << _weight << endl;
-       return _weight;
-}
-
-
-// calculates complete event from the effective masses of the (i + 1)-body
-// systems, the Lorentz vector of the decaying system, and the decay angles
-// uses the break-up momenta calculated by calcWeight()
-void
-nBodyPhaseSpaceGen::calcEventKinematics(const lorentzVector& nBody)  // Lorentz vector of n-body system in lab frame
-{
-       // build event starting in n-body RF going down to 2-body RF
-       // is more efficicient than Raubold-Lynch method, since it requitres only one rotation and boost per daughter
-       lorentzVector P = nBody;  // Lorentz of (i + 1)-body system in lab frame
-       for (unsigned int i = _n - 1; i >= 1; --i) {  // loop from n-body down to 2-body
-               // construct Lorentz vector of daughter _m[i] in (i + 1)-body RF
-               const double   sinTheta = sqrt(1 - _cosTheta[i] * _cosTheta[i]);
-               const double   pT       = _breakupMom[i] * sinTheta;
-               lorentzVector& daughter = _daughters[i];
-               daughter.SetPxPyPzE(pT * cos(_phi[i]),
-                                   pT * sin(_phi[i]),
-                                   _breakupMom[i] * _cosTheta[i],
-                                   sqrt(_m[i] * _m[i] + _breakupMom[i] * _breakupMom[i]));
-               // boost daughter into lab frame
-               daughter.Boost(P.BoostVector());
-               // calculate Lorentz vector of i-body system in lab frame
-               P -= daughter;
-       }
-       // set last daughter
-       _daughters[0] = P;
-}
-
-
-// calculates maximum weight for given n-body mass
-double
-nBodyPhaseSpaceGen::estimateMaxWeight(const double       nBodyMass,        // sic!
-                                      const unsigned int nmbOfIterations)  // number of generated events
-{
-       double maxWeight = 0;
-       for (unsigned int i = 0; i < nmbOfIterations; ++i) {
-               pickMasses(nBodyMass);
-               calcWeight();
-               maxWeight = max(_weight, maxWeight);
-       }
-       return maxWeight;
-}
-
-
-ostream&
-nBodyPhaseSpaceGen::print(ostream& out) const
-{
-       out << "nBodyPhaseSpaceGen parameters:" << endl
-           << "    number of daughter particles ............... " << _n                 << endl
-           << "    masses of the daughter particles ........... " << _m                 << endl
-           << "    sums of daughter particle masses ........... " << _mSum              << endl
-           << "    effective masses of (i + 1)-body systems ... " << _M                 << endl
-           << "    cos(polar angle) in (i + 1)-body systems ... " << _cosTheta          << endl
-           << "    azimuth in (i + 1)-body systems ............ " << _phi               << endl
-           << "    breakup momenta in (i + 1)-body systems .... " << _breakupMom        << endl
-           << "    normalization value ........................ " << _norm              << endl
-           << "    weight of generated event .................. " << _weight            << endl
-           << "    maximum weight used in hit-miss MC ......... " << _maxWeight         << endl
-           << "    maximum weight since instantiation ......... " << _maxWeightObserved << endl
-           << "    daughter four-momenta:" << endl;
-       for (unsigned int i = 0; i < _n; ++i)
-               out << "        daughter " << i << ": " << _daughters[i] << endl;
-       return out;
-}
diff --git a/STARLIGHT/starlight/src/nucleus.cpp~ b/STARLIGHT/starlight/src/nucleus.cpp~
deleted file mode 100644 (file)
index d175246..0000000
+++ /dev/null
@@ -1,216 +0,0 @@
-///////////////////////////////////////////////////////////////////////////
-//
-//    Copyright 2010
-//
-//    This file is part of starlight.
-//
-//    starlight 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 3 of the License, or
-//    (at your option) any later version.
-//
-//    starlight 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 starlight. If not, see <http://www.gnu.org/licenses/>.
-//
-///////////////////////////////////////////////////////////////////////////
-//
-// File and Version Information:
-// $Rev:: 164                         $: revision of last commit
-// $Author:: odjuvsla                 $: author of last commit
-// $Date:: 2013-10-06 16:18:08 +0200 #$: date of last commit
-//
-// Description:
-//
-//
-//
-///////////////////////////////////////////////////////////////////////////
-
-
-#include <iostream>
-#include <fstream>
-
-#include "starlightconstants.h"
-#include "reportingUtils.h"
-#include "nucleus.h"
-#include <inputParameters.h>
-
-
-using namespace std;
-using namespace starlightConstants;
-
-//______________________________________________________________________________
-nucleus::nucleus(const int    Z,
-                 const int    A,
-                 const double deuteronSlopePar,
-                 const bool   dAuCoherentProduction)
-       : _Z(Z),
-         _A(A),
-         _deuteronSlopePar(deuteronSlopePar),
-         _dAuCoherentProduction(dAuCoherentProduction)
-{
-  init();      
-}
-
-void nucleus::init()
-{
-  switch (_Z) {
-       case 79:
-               {
-                       _Q0   = 0.060;
-                       _rho0 = 0.159407;
-               }
-               break;
-       case 53:
-       case 49:
-               {
-                       _Q0   = 0.069;
-                       _rho0 = 0.161626;
-               }
-               break;
-       case 29:
-               {
-                       _Q0   = 0.087;
-                       _rho0 = 0.166878;
-               }
-               break;
-       case 14:
-               {
-                       _Q0   = 0.115;
-                       _rho0 = 0.177128;
-               }
-               break;
-       case 8:
-               {
-                       _Q0   = 0.138;
-                       _rho0 = 0.188459;
-               }
-               break;
-       case 82:
-               {
-                       _Q0   = 0.059;
-                       _rho0 = 0.159176;
-               }
-               break;
-       case 20:
-               {
-                       _Q0   = 0.102;
-                       _rho0 = 0.171907;
-               }
-               break;
-       case 1: // _Q0 and _rho0 are not relevant for protons.
-               {
-                       _Q0   = -1.0;
-                       _rho0 = -1.0;
-               }
-               break;
-       default:
-               printWarn << "density not defined for projectile with Z = " << _Z << ". using defaults." << endl;
-               _rho0 = 0.16;  //'typical' density
-               _Q0   = 0.0;   // not used.
-       }
-       _r0 = 1.16 * (1. - 1.16 * pow(_A, -2. / 3.));  // for FRITIOF and FormFactor.
-}
-
-//______________________________________________________________________________
-nucleus::~nucleus()
-{ }
-
-
-//______________________________________________________________________________
-double
-nucleus::nuclearRadius() const
-{
-       // we use this for specific nuclei, Au, Pb, protons...
-       if (_Z == 79)                // Au
-               return 6.38;  // [fm]
-       if (_Z == 82)                // Pb
-               return 6.62;  // [fm]
-       if ((_Z == 1) && (_A == 1))  // proton
-               return 0.7;   // [fm]
-       if ((_Z == 1) && (_A == 2))  // deuteron
-               return 1.9;   // [fm]
-       return 1.2 * pow(_A, 1. / 3.); 
-}
-
-
-//______________________________________________________________________________
-double
-nucleus::formFactor(const double t) const
-{
-       // electromagnetic form factor of proton
-       if ((_Z == 1) && (_A == 1)) {
-               const double rec = 1. / (1. + t / 0.71);
-               return rec * rec;
-       }
-       // deuteron form factor
-       if ((_Z == 1) && (_A == 2)) {   // careful with this line on dAu
-               // this is for dAu//Sergey
-               // sergey's stuff, also replaced b with _deuteronSlopePar and dropped t02 since it wasnt used
-               // incoherent form factor F(t) = 0.34 e(141.5 t) + 0.58 e(26.1 t) + 0.08 e(15.5 t)
-               const double st  = 0.34 * exp(-141.5 * t    ) + 0.58 * exp(-26.1 * t    ) + 0.08 * exp(-15.5 * t    );
-               const double st4 = 0.34 * exp(-141.5 * t / 4) + 0.58 * exp(-26.1 * t / 4) + 0.08 * exp(-15.5 * t / 4);
-               // st paramters from Franco and Varma for st eqn PRL33 ...
-               // form factor from Eisenberg, nuclear physics B 104
-               const double arg = _deuteronSlopePar * t;
-               if (_dAuCoherentProduction)
-                       return (st4 * st4 * exp(-arg) - 0.068 * st4 * exp(-arg * 3. / 4.));
-               return exp(-arg) * 0.5 * (1 + st) - 0.068 * exp(-arg * 3. / 4.)
-                       - st4 * st4 * exp(-arg) + 0.068 * st4 * exp(-arg * 3. / 4.);
-       }
-       // nuclear form factor
-       // use parameterization from FRITIOF
-       // R = r0 * A^{1/3} with r0 = 1.16 * (1 - 1.16 * A^{-2/3})
-       const double R    = fritiofR0();
-       const double q    = sqrt(t);
-       const double arg1 = q * R / hbarc;
-       const double arg2 = hbarc / (q * _r0);
-       const double sph  = (sin(arg1) - arg1 * cos(arg1)) * 3. * arg2 * arg2 * arg2 / double(_A);
-       const double a0   = 0.70;  // [fm]
-       return sph / (1. + (a0 * a0 * t) / (hbarc * hbarc));
-}
-
-//______________________________________________________________________________
-
-double
-nucleus::dipoleFormFactor(const double t, const double t0) const
-{
-     const double rec = 1. / (1. + t / t0);
-     return rec * rec;
-}
-
-//______________________________________________________________________________
-double
-nucleus::thickness(const double b) const
-{
-       //    JS      This code calculates the nuclear thickness function as per Eq. 4 in
-       //    Klein and Nystrand, PRC 60.
-       //    former DOUBLE PRECISION FUNCTION T(b)
-                                                  
-       // data for Gauss integration
-       const unsigned int nmbPoints         = 5;
-       const double       xg[nmbPoints + 1] = {0., 0.1488743390, 0.4333953941, 0.6794095683,
-                                               0.8650633667, 0.9739065285};
-       const double       ag[nmbPoints + 1] = {0., 0.2955242247, 0.2692667193, 0.2190863625,
-                                               0.1494513492, 0.0666713443};
-  
-       const double zMin   = 0;
-       const double zMax   = 15;
-       const double zRange = 0.5 * (zMax - zMin); 
-       const double zMean  = 0.5 * (zMax + zMin); 
-       double       sum    = 0;
-       for(unsigned int i = 1; i <= nmbPoints; ++i) {
-               double zsp    = zRange * xg[i] + zMean;
-               double radius = sqrt(b * b + zsp * zsp);
-               sum          += ag[i] * rws(radius);
-               zsp           = zRange * (-xg[i]) + zMean;
-               radius        = sqrt(b * b + zsp * zsp);
-               sum          += ag[i] * rws(radius);
-       }
-  
-       return 2. * zRange * sum;
-}
diff --git a/STARLIGHT/starlight/src/wideResonanceCrossSection.cpp~ b/STARLIGHT/starlight/src/wideResonanceCrossSection.cpp~
deleted file mode 100644 (file)
index 3a43ffa..0000000
+++ /dev/null
@@ -1,171 +0,0 @@
-///////////////////////////////////////////////////////////////////////////
-//
-//    Copyright 2010
-//
-//    This file is part of starlight.
-//
-//    starlight 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 3 of the License, or
-//    (at your option) any later version.
-//
-//    starlight 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 starlight. If not, see <http://www.gnu.org/licenses/>.
-//
-///////////////////////////////////////////////////////////////////////////
-//
-// File and Version Information:
-// $Rev:: 164                         $: revision of last commit
-// $Author:: odjuvsla                 $: author of last commit
-// $Date:: 2013-10-06 16:18:08 +0200 #$: date of last commit
-//
-// Description:
-//
-//
-//
-///////////////////////////////////////////////////////////////////////////
-
-
-#include <iostream>
-#include <fstream>
-#include <cmath>
-
-#include "starlightconstants.h"
-#include "wideResonanceCrossSection.h"
-
-
-using namespace std;
-using namespace starlightConstants;
-
-
-//______________________________________________________________________________
-wideResonanceCrossSection::wideResonanceCrossSection(const beamBeamSystem&  bbsystem)
-       : photonNucleusCrossSection(bbsystem)//hrm
-{
-       _wideWmax = _wMax;
-       _wideWmin = _wMin;
-       _wideYmax = _yMax;
-       _wideYmin = -1.0 * _wideYmax;
-       _Ep       = inputParametersInstance.protonEnergy();
-}
-
-
-//______________________________________________________________________________
-wideResonanceCrossSection::~wideResonanceCrossSection()
-{
-
-}
-
-
-//______________________________________________________________________________
-void
-wideResonanceCrossSection::crossSectionCalculation(const double bwnormsave)
-{
-       //     This subroutine calculates the cross-section assuming a wide
-       //     (Breit-Wigner) resonance.
-
-       // double Av,Wgp,cs,cvma;
-       double W,dW,dY;
-       double y1,y2,y12,ega1,ega2,ega12;
-       // double t,tmin,tmax;
-       double csgA1,csgA2,csgA12,int_r,dR,rate;
-       double dsigdW,dsigdWalt,dndW,tmp;
-       double dsigdW2;
-       // double ax,bx;
-       double Eth;
-       int    I,J,NW,NY;
-       // int    K,NGAUSS;
-                                                                                                                                                      
-       // ----------------- !!!!!!!!!!!!!!!!!!!! -----------------------------
-                                                                                                                                                      
-       double bwnorm =bwnormsave;//used to transfer the bwnorm from the luminosity tables
-
-       // --------------------------------------------------------------------
-       //gamma+nucleon threshold.
-
-       Eth=0.5*(((_wideWmin+protonMass)*(_wideWmin+protonMass)
-                 -protonMass*protonMass)/(_Ep+sqrt(_Ep*_Ep-protonMass*protonMass)));
-                                                                                                                                                      
-       NW   = 100;
-       dW   = (_wideWmax-_wideWmin)/double(NW);
-  
-       NY   =  1200;
-       dY   = (_wideYmax-_wideYmin)/double(NY);
-  
-       if (getBNORM()  ==  0.){
-               cout<<" Using Breit-Wigner Resonance Profile."<<endl;
-       }
-       else{
-               cout<<" Using Breit-Wigner plus direct pi+pi- profile."<<endl;
-       }
-  
-       cout<<" Integrating over W from "<<_wideWmin<<" to "<<_wideWmax<<endl;
-                                                                                                                                                      
-       int_r=0.;
-       for(I=0;I<=NW-1;I++){
-    
-               W = _wideWmin + double(I)*dW + 0.5*dW;
-    
-               tmp = 0.0;
-               dsigdW=0.0;
-               dsigdW2=0.0;
-               dsigdWalt=0.0;
-               dndW=0.0;
-    
-               for(J=0;J<=NY-1;J++){
-      
-                       y1  = _wideYmin + double(J)*dY;
-                       y2  = _wideYmin + double(J+1)*dY;
-                       y12 = 0.5*(y1+y2);
-      
-                       ega1  = 0.5*W*exp(y1);
-                       ega2  = 0.5*W*exp(y2);
-                       ega12 = 0.5*W*exp(y12);
-      
-                       if(ega1 < Eth) continue;
-                       if(ega2 > maxPhotonEnergy()) continue;
-                       // check it !!
-          
-                       if(J == 0){
-                               // >> 1st Point (Calculated only the first time)     =====>>>
-                               //ega1 used.                                                        
-                               csgA1=getcsgA(ega1,W);
-                       }
-                       else{
-                               csgA1 = csgA2;
-                       }
-          
-                       //         >> Middle Point                      =====>>>
-                       csgA12=getcsgA(ega12,W);         
-
-                       //         >> Second Point                      =====>>>
-                       csgA2=getcsgA(ega2,W);
-      
-                       //>> Sum the contribution for this W,Y. The 2 accounts for the 2 beams
-                       dR  = ega1*photonFlux(ega1)*csgA1;
-                       dR  = dR + 4.*ega12*photonFlux(ega12)*csgA12;
-                       dR  = dR + ega2*photonFlux(ega2)*csgA2;
-                       tmp = tmp+2.*dR*(dY/6.);
-                       dR  = dR*(dY/6.)*breitWigner(W,bwnorm)*dW;
-      
-                       //For identical beams, we double.  Either may emit photon/scatter
-                       //For large differences in Z, we approx, that only beam1 emits photon
-                       //and beam2 scatters, eg d-Au where beam1=au and beam2=d
-                       if(getbbs().beam1().A()==getbbs().beam2().A()){
-                               dR  = 2.*dR;
-                       }
-                       int_r = int_r+dR;  
-               }
-       }
-                                                                                                                                                      
-       rate=luminosity()*int_r;
-  
-       cout<<" Cross section (mb): "<<10.*int_r<<endl;
-       cout<<" Production rate   : "<<rate<<" Hz"<<endl;
-}
diff --git a/STARLIGHT/test/Config.C~ b/STARLIGHT/test/Config.C~
deleted file mode 100644 (file)
index 7fa5431..0000000
+++ /dev/null
@@ -1,472 +0,0 @@
-//
-// Configuration for p-p Diffraction studies
-// STARLIGHT pp 7000
-
-// One can use the configuration macro in compiled mode by
-// root [0] gSystem->Load("libgeant321");
-// root [0] gSystem->SetIncludePath("-I$ROOTSYS/include -I$ALICE_ROOT/include\
-//                   -I$ALICE_ROOT -I$ALICE/geant3/TGeant3");
-// root [0] .x grun.C(1,"Config.C++")
-
-#if !defined(__CINT__) || defined(__MAKECINT__)
-#include <Riostream.h>
-#include <TRandom.h>
-#include <TDatime.h>
-#include <TSystem.h>
-#include <TVirtualMC.h>
-#include <TGeant3TGeo.h>
-#include "STEER/AliRunLoader.h"
-#include "STEER/AliRun.h"
-#include "STEER/AliConfig.h"
-#include "PYTHIA6/AliDecayerPythia.h"
-#include "PYTHIA6/AliGenPythia.h"
-#include "TDPMjet/AliGenDPMjet.h"
-#include "STEER/AliMagFCheb.h"
-#include "STRUCT/AliBODY.h"
-#include "STRUCT/AliMAG.h"
-#include "STRUCT/AliABSOv3.h"
-#include "STRUCT/AliDIPOv3.h"
-#include "STRUCT/AliHALLv3.h"
-#include "STRUCT/AliFRAMEv2.h"
-#include "STRUCT/AliSHILv3.h"
-#include "STRUCT/AliPIPEv3.h"
-#include "ITS/AliITSv11Hybrid.h"
-#include "TPC/AliTPCv2.h"
-#include "TOF/AliTOFv6T0.h"
-#include "HMPID/AliHMPIDv3.h"
-#include "ZDC/AliZDCv3.h"
-#include "TRD/AliTRDv1.h"
-#include "TRD/AliTRDgeometry.h"
-#include "FMD/AliFMDv1.h"
-#include "MUON/AliMUONv1.h"
-#include "PHOS/AliPHOSv1.h"
-#include "PHOS/AliPHOSSimParam.h"
-#include "PMD/AliPMDv1.h"
-#include "T0/AliT0v1.h"
-#include "EMCAL/AliEMCALv2.h"
-#include "ACORDE/AliACORDEv1.h"
-#include "VZERO/AliVZEROv7.h"
-#endif
-
-
-enum PDC06Proc_t
-{
-  kStarlight, kRunMax,
-};
-
-const char * pprRunName[] = {
-  "kStarlight"
-};
-
-enum Mag_t
-{
-  kNoField, k5kG, kFieldMax
-};
-
-const char * pprField[] = {
-  "kNoField", "k5kG"
-};
-
-//--- Functions ---
-void ProcessEnvironmentVars();
-
-// Generator, field, beam energy
-static PDC06Proc_t   proc         = kStarlight;
-static Mag_t         mag          = k5kG;
-static Float_t       energy       = 7000; // energy in CMS
-static Float_t       pBeamEnergy  = 3500.; // energy p-Beam
-static Int_t         runNumber    = 0; 
-static Int_t         eventOffset  = 0;
-
-//========================//
-// Set Random Number seed //
-//========================//
-TDatime dt;
-static UInt_t seed    = dt.Get();
-
-// Comment line
-static TString comment;
-
-void Config()
-{
-  // Get settings from environment variables
-  ProcessEnvironmentVars();
-
-  gRandom->SetSeed(seed);
-  cerr<<"Seed for random number generation= "<<seed<<endl;
-
-  // Libraries required by geant321
-#if defined(__CINT__)
-  gSystem->Load("liblhapdf");      // Parton density functions
-  gSystem->Load("libEGPythia6");   // TGenerator interface
-//   if (proc == kPythia6 || proc == kPhojet) {
-//     gSystem->Load("libpythia6");        // Pythia 6.2
-//   } else {
-     gSystem->Load("libpythia6.4.21");   // Pythia 6.4
-//   }
-  gSystem->Load("libAliPythia6");  // ALICE specific implementations
-  gSystem->Load("libgeant321");
-#endif
-
-  new TGeant3TGeo("C++ Interface to Geant3");
-
-  //=======================================================================
-  //  Create the output file
-
-
-  AliRunLoader* rl=0x0;
-
-  cout<<"Config.C: Creating Run Loader ..."<<endl;
-  rl = AliRunLoader::Open("galice.root",
-                         AliConfig::GetDefaultEventFolderName(),
-                         "recreate");
-  if (rl == 0x0)
-    {
-      gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
-      return;
-    }
-  rl->SetCompressionLevel(2);
-  rl->SetNumberOfEventsPerFile(1000);
-  gAlice->SetRunLoader(rl);
-  // gAlice->SetGeometryFromFile("geometry.root");
-  // gAlice->SetGeometryFromCDB();
-
-  // Set the trigger configuration: proton-proton
-  AliSimulation::Instance()->SetTriggerConfig("p-p"); 
-
-  //
-  //=======================================================================
-  // ************* STEERING parameters FOR ALICE SIMULATION **************
-  // --- Specify event type to be tracked through the ALICE setup
-  // --- All positions are in cm, angles in degrees, and P and E in GeV
-
-  gMC->SetProcess("DCAY",1);
-  gMC->SetProcess("PAIR",1);
-  gMC->SetProcess("COMP",1);
-  gMC->SetProcess("PHOT",1);
-  gMC->SetProcess("PFIS",0);
-  gMC->SetProcess("DRAY",0);
-  gMC->SetProcess("ANNI",1);
-  gMC->SetProcess("BREM",1);
-  gMC->SetProcess("MUNU",1);
-  gMC->SetProcess("CKOV",1);
-  gMC->SetProcess("HADR",1);
-  gMC->SetProcess("LOSS",2);
-  gMC->SetProcess("MULS",1);
-  gMC->SetProcess("RAYL",1);
-
-  Float_t cut = 1.e-3;        // 1MeV cut by default
-  Float_t tofmax = 1.e10;
-
-  gMC->SetCut("CUTGAM", cut);
-  gMC->SetCut("CUTELE", cut);
-  gMC->SetCut("CUTNEU", cut);
-  gMC->SetCut("CUTHAD", cut);
-  gMC->SetCut("CUTMUO", cut);
-  gMC->SetCut("BCUTE",  cut);
-  gMC->SetCut("BCUTM",  cut);
-  gMC->SetCut("DCUTE",  cut);
-  gMC->SetCut("DCUTM",  cut);
-  gMC->SetCut("PPCUTM", cut);
-  gMC->SetCut("TOFMAX", tofmax);
-
-  //======================//
-  // Set External decayer //
-  //======================//
-  TVirtualMCDecayer* decayer = new AliDecayerPythia();
-  decayer->SetForceDecay(kAll);
-
-  decayer->Init();
-  gMC->SetExternalDecayer(decayer);
-
-  //=========================//
-  // Generator Configuration //
-  //=========================//
-  AliGenerator* gener = 0x0;
-
-  if (proc == kStarlight) {
-    gSystem->SetIncludePath("-I./ -I$ROOTSYS/include -I$ALICE_ROOT/include -I$ALICE_ROOT/EVGEN -I$ALICE_ROOT -I$ALICE/geant3/TGeant3");
-    gROOT->LoadMacro("AliGenExtFile_.cxx+");
-    gROOT->LoadMacro("AliGenReaderSL_.cxx+");
-
-    AliGenExtFile_* genExtFile = new AliGenExtFile_(-1);
-    AliGenReaderSL_* reader = new AliGenReaderSL_;
-    reader->SetFileName("pp_to_f0.out");
-    reader->SetFormat(3);
-    genExtFile->SetReader(reader);
-    genExtFile->SetStartEvent(eventOffset); // use events starting from eventOffset
-    gener= genExtFile;
-  }
-
-  //
-  //
-  // Size of the interaction diamond
-  // Longitudinal
-  Float_t sigmaz  = 5.4 / TMath::Sqrt(2.); // [cm]
-  //   if (energy == 900)
-  //sigmaz  = 10.5 / TMath::Sqrt(2.); // [cm]
-  //sigmaz = 3.7;
-  if (energy == 7000)
-    sigmaz  = 6.3 / TMath::Sqrt(2.); // [cm]
-  
-  //
-  // Transverse
-  Float_t betast  = 10;                 // beta* [m]
-  if (runNumber >= 117048) betast = 2.;
-  if (runNumber >  122375) betast = 3.5; // starting with fill 1179
-
-  Float_t eps     = 5.e-6;            // emittance [m]
-  Float_t gamma   = energy / 2.0 / 0.938272;  // relativistic gamma [1]
-  Float_t sigmaxy = TMath::Sqrt(eps * betast / gamma) / TMath::Sqrt(2.) * 100.;  // [cm]
-  printf("\n \n Diamond size x-y: %10.3e z: %10.3e\n \n", sigmaxy, sigmaz);
-
-  gener->SetSigma(sigmaxy, sigmaxy, sigmaz);      // Sigma in (X,Y,Z) (cm) on IP position
-  gener->SetVertexSmear(kPerEvent);
-  gener->Init();
-
-  printf("\n \n Comment: %s \n \n", comment.Data());
-
-  //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-  TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., AliMagF::k5kG,
-                                                      AliMagF::kBeamTypepp, pBeamEnergy));
-
-  rl->CdGAFile();
-
-  Int_t iABSO  = 1;
-  Int_t iACORDE= 0;
-  Int_t iDIPO  = 1;
-  Int_t iEMCAL = 1;
-  Int_t iFMD   = 1;
-  Int_t iFRAME = 1;
-  Int_t iHALL  = 1;
-  Int_t iITS   = 1;
-  Int_t iMAG   = 1;
-  Int_t iMUON  = 1;
-  Int_t iPHOS  = 1;
-  Int_t iPIPE  = 1;
-  Int_t iPMD   = 1;
-  Int_t iHMPID = 1;
-  Int_t iSHIL  = 1;
-  Int_t iT0    = 1;
-  Int_t iTOF   = 1;
-  Int_t iTPC   = 1;
-  Int_t iTRD   = 1;
-  Int_t iVZERO = 1;
-  Int_t iZDC   = 1;
-
-
-  //=================== Alice BODY parameters =============================
-  AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
-
-
-  if (iMAG)
-    {
-      //=================== MAG parameters ============================
-      // --- Start with Magnet since detector layouts may be depending ---
-      // --- on the selected Magnet dimensions ---
-      AliMAG *MAG = new AliMAG("MAG", "Magnet");
-    }
-
-
-  if (iABSO)
-    {
-      //=================== ABSO parameters ============================
-      AliABSO *ABSO = new AliABSOv3("ABSO", "Muon Absorber");
-    }
-
-  if (iDIPO)
-    {
-      //=================== DIPO parameters ============================
-
-      AliDIPO *DIPO = new AliDIPOv3("DIPO", "Dipole version 3");
-    }
-
-  if (iHALL)
-    {
-      //=================== HALL parameters ============================
-
-      AliHALL *HALL = new AliHALLv3("HALL", "Alice Hall");
-    }
-
-
-  if (iFRAME)
-    {
-      //=================== FRAME parameters ============================
-
-      AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
-      FRAME->SetHoles(1);
-    }
-
-  if (iSHIL)
-    {
-      //=================== SHIL parameters ============================
-
-      AliSHIL *SHIL = new AliSHILv3("SHIL", "Shielding Version 3");
-    }
-
-
-  if (iPIPE)
-    {
-      //=================== PIPE parameters ============================
-
-      AliPIPE *PIPE = new AliPIPEv3("PIPE", "Beam Pipe");
-    }
-
-  if (iITS)
-    {
-      //=================== ITS parameters ============================
-
-      AliITS *ITS  = new AliITSv11("ITS","ITS v11");
-    }
-
-  if (iTPC)
-    {
-      //============================ TPC parameters =====================
-
-      AliTPC *TPC = new AliTPCv2("TPC", "Default");
-    }
-
-
-  if (iTOF) {
-    //=================== TOF parameters ============================
-
-    AliTOF *TOF = new AliTOFv6T0("TOF", "normal TOF");
-  }
-
-
-  if (iHMPID)
-    {
-      //=================== HMPID parameters ===========================
-
-      AliHMPID *HMPID = new AliHMPIDv3("HMPID", "normal HMPID");
-
-    }
-
-
-  if (iZDC)
-    {
-      //=================== ZDC parameters ============================
-
-      AliZDC *ZDC = new AliZDCv3("ZDC", "normal ZDC");
-    }
-
-  if (iTRD)
-    {
-      //=================== TRD parameters ============================
-
-      AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
-      AliTRDgeometry *geoTRD = TRD->GetGeometry();
-      // Partial geometry: modules at 0,1,7,8,9,16,17
-      // starting at 3h in positive direction
-      geoTRD->SetSMstatus(2,0);
-      geoTRD->SetSMstatus(3,0);
-      geoTRD->SetSMstatus(4,0);
-      geoTRD->SetSMstatus(5,0);
-      geoTRD->SetSMstatus(6,0);
-      geoTRD->SetSMstatus(11,0);
-      geoTRD->SetSMstatus(12,0);
-      geoTRD->SetSMstatus(13,0);
-      geoTRD->SetSMstatus(14,0);
-      geoTRD->SetSMstatus(15,0);
-      geoTRD->SetSMstatus(16,0);
-    }
-
-  if (iFMD)
-    {
-      //=================== FMD parameters ============================
-
-      AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
-    }
-
-  if (iMUON)
-    {
-      //=================== MUON parameters ===========================
-      // New MUONv1 version (geometry defined via builders)
-
-      AliMUON *MUON = new AliMUONv1("MUON", "default");
-    }
-
-  if (iPHOS)
-    {
-      //=================== PHOS parameters ===========================
-
-      AliPHOS *PHOS = new AliPHOSv1("PHOS", "noCPV_Modules123");
-
-    }
-
-
-  if (iPMD)
-    {
-      //=================== PMD parameters ============================
-
-      AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
-    }
-
-  if (iT0)
-    {
-      //=================== T0 parameters ============================
-      AliT0 *T0 = new AliT0v1("T0", "T0 Detector");
-    }
-
-  if (iEMCAL)
-    {
-      //=================== EMCAL parameters ============================
-
-      AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "EMCAL_FIRSTYEAR");
-    }
-
-  if (iACORDE)
-    {
-      //=================== ACORDE parameters ============================
-
-      AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE");
-    }
-
-  if (iVZERO)
-    {
-      //=================== ACORDE parameters ============================
-
-      AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO");
-    }
-}
-
-void ProcessEnvironmentVars()
-{
-  // Run type
-  if (gSystem->Getenv("CONFIG_RUN_TYPE")) {
-    for (Int_t iRun = 0; iRun < kRunMax; iRun++) {
-      if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun])==0) {
-       proc = (PDC06Proc_t)iRun;
-       cout<<"Run type set to "<<pprRunName[iRun]<<endl;
-      }
-    }
-  }
-  // Field
-  if (gSystem->Getenv("CONFIG_FIELD")) {
-    for (Int_t iField = 0; iField < kFieldMax; iField++) {
-      if (strcmp(gSystem->Getenv("CONFIG_FIELD"), pprField[iField])==0) {
-       mag = (Mag_t)iField;
-       cout<<"Field set to "<<pprField[iField]<<endl;
-      }
-    }
-  }
-  // Energy
-  if (gSystem->Getenv("CONFIG_ENERGY")) {
-    energy = atoi(gSystem->Getenv("CONFIG_ENERGY"));
-    cout<<"Energy set to "<<energy<<" GeV"<<endl;
-  }
-  // Random Number seed
-  if (gSystem->Getenv("CONFIG_SEED")) {
-    seed = atoi(gSystem->Getenv("CONFIG_SEED"));
-  }
-  // Run number
-  if (gSystem->Getenv("DC_RUN")) {
-    runNumber = atoi(gSystem->Getenv("DC_RUN"));
-  }
-  // Event offset
-  if (gSystem->Getenv("DC_EVENT")) {
-    eventOffset = 400*(atoi(gSystem->Getenv("DC_EVENT")) - 1); // 400 events per job (->sim.C)
-    std::cout << "eventOffset= " << eventOffset << std::endl;
-  }
-}
-
diff --git a/STARLIGHT/test/Config_.C~ b/STARLIGHT/test/Config_.C~
deleted file mode 100644 (file)
index fdb3e02..0000000
+++ /dev/null
@@ -1,384 +0,0 @@
-//Configuration of simulation
-
-enum PprRun_t 
-{
-    kSTARLIGHT,
-    kRunMax
-};
-
-const char* pprRunName[] = {
-  "kSTARLIGHT"
-};
-
-enum PprTrigConf_t
-{
-    kDefaultPPTrig, kDefaultPbPbTrig
-};
-
-const char * pprTrigConfName[] = {
-    "p-p","Pb-Pb"
-};
-
-// This part for configuration    
-
-static PprRun_t srun = kSTARLIGHT;
-static AliMagF::BMap_t smag = AliMagF::k5kG;
-static Int_t    sseed = 12345; //Set 0 to use the current time
-static PprTrigConf_t strig = kDefaultPPTrig; // default pp trigger configuration
-
-// Comment line 
-static TString  comment;
-
-// Functions
-Float_t EtaToTheta(Float_t arg);
-void ProcessEnvironmentVars();
-
-void Config()
-{
-    // ThetaRange is (0., 180.). It was (0.28,179.72) 7/12/00 09:00
-    // Theta range given through pseudorapidity limits 22/6/2001
-
-    // Get settings from environment variables
-    ProcessEnvironmentVars();
-
-    // Set Random Number seed
-    gRandom->SetSeed(sseed);
-    cout<<"Seed for random number generation= "<<gRandom->GetSeed()<<endl; 
-
-
-   // libraries required by geant321 and Pythia: loaded in sim.C
-
-    new     TGeant3TGeo("C++ Interface to Geant3");
-
-  // Output every 100 tracks
-
-    TVirtualMC * vmc = TVirtualMC::GetMC();
-
-//   ((TGeant3*)vmc)->SetSWIT(4,100);
-
-    AliRunLoader* rl=0x0;
-
-    AliLog::Message(AliLog::kInfo, "Creating Run Loader", "", "", "Config()"," ConfigPPR.C", __LINE__);
-
-    rl = AliRunLoader::Open("galice.root",
-                           AliConfig::GetDefaultEventFolderName(),
-                           "recreate");
-    if (rl == 0x0)
-      {
-       gAlice->Fatal("Config.C","Can not instatiate the Run Loader");
-       return;
-      }
-    rl->SetCompressionLevel(2);
-    rl->SetNumberOfEventsPerFile(100);
-    gAlice->SetRunLoader(rl);
-
-    // Set the trigger configuration
-    AliSimulation::Instance()->SetTriggerConfig(pprTrigConfName[strig]);
-    cout<<"Trigger configuration is set to  "<<pprTrigConfName[strig]<<endl;
-
-    //
-    // Set External decayer
-    AliDecayer *decayer = new AliDecayerPythia();
-
-    decayer->Init();
-    vmc->SetExternalDecayer(decayer);
-    //
-    //
-    //=======================================================================
-    //
-    //=======================================================================
-    // ************* STEERING parameters FOR ALICE SIMULATION **************
-    // --- Specify event type to be tracked through the ALICE setup
-    // --- All positions are in cm, angles in degrees, and P and E in GeV
-
-    vmc->SetProcess("DCAY",1);
-    vmc->SetProcess("PAIR",1);
-    vmc->SetProcess("COMP",1);
-    vmc->SetProcess("PHOT",1);
-    vmc->SetProcess("PFIS",0);
-    vmc->SetProcess("DRAY",0);
-    vmc->SetProcess("ANNI",1);
-    vmc->SetProcess("BREM",1);
-    vmc->SetProcess("MUNU",1);
-    vmc->SetProcess("CKOV",1);
-    vmc->SetProcess("HADR",1);
-    vmc->SetProcess("LOSS",2);
-    vmc->SetProcess("MULS",1);
-    vmc->SetProcess("RAYL",1);
-
-    Float_t cut = 1.e-3;        // 1MeV cut by default
-    Float_t tofmax = 1.e10;
-
-    vmc->SetCut("CUTGAM", cut);
-    vmc->SetCut("CUTELE", cut);
-    vmc->SetCut("CUTNEU", cut);
-    vmc->SetCut("CUTHAD", cut);
-    vmc->SetCut("CUTMUO", cut);
-    vmc->SetCut("BCUTE",  cut); 
-    vmc->SetCut("BCUTM",  cut); 
-    vmc->SetCut("DCUTE",  cut); 
-    vmc->SetCut("DCUTM",  cut); 
-    vmc->SetCut("PPCUTM", cut);
-    vmc->SetCut("TOFMAX", tofmax); 
-
-    // Generator Configuration
-    gSystem->AddDynamicPath("../../objdir/lib/tgt_linux");
-    gSystem->Load("libStarLight");
-    gSystem->Load("libAliStarLight.so");
-
-    AliGenStarLight* sl = new AliGenStarLight(1000*1000);
-
-    sl->SetParameter("BEAM_1_Z     =   82    #Z of projectile");
-    sl->SetParameter("BEAM_1_A     =  208    #A of projectile");
-    sl->SetParameter("BEAM_2_Z     =   82    #Z of target");
-    sl->SetParameter("BEAM_2_A     =  208    #A of target");
-    sl->SetParameter("BEAM_1_GAMMA = 1470    #Gamma of the colliding ions");
-    sl->SetParameter("BEAM_2_GAMMA = 1470    #Gamma of the colliding ions");
-    sl->SetParameter("W_MAX        =   12.0  #Max value of w");
-    sl->SetParameter("W_MIN        =    2.0  #Min value of w");
-    sl->SetParameter("W_N_BINS     =   40    #Bins i w");
-    sl->SetParameter("RAP_MAX      =    8.   #max y");
-    sl->SetParameter("RAP_N_BINS   =   80    #Bins i y");
-    sl->SetParameter("CUT_PT       =    0    #Cut in pT? 0 = (no, 1 = yes)");
-    sl->SetParameter("PT_MIN       =    1.0  #Minimum pT in GeV");
-    sl->SetParameter("PT_MAX       =    3.0  #Maximum pT in GeV");
-    sl->SetParameter("CUT_ETA      =    0    #Cut in pseudorapidity? (0 = no, 1 = yes)");
-    sl->SetParameter("ETA_MIN      =  -10    #Minimum pseudorapidity");
-    sl->SetParameter("ETA_MAX      =   10    #Maximum pseudorapidity");
-    sl->SetParameter("PROD_MODE    =    2    #gg or gP switch (1 = 2-photon, 2 = coherent vector meson (narrow), 3 = coherent vector meson (wide), # 4 = incoherent vector meson, 5 = A+A DPMJet single, 6 = A+A DPMJet double, 7 = p+A DPMJet single, 8 = p+A Pythia single )");
-    // is N_EVENTS valid
-    sl->SetParameter("N_EVENTS     = 10000000  #Number of events");
-    sl->SetParameter("PROD_PID     =  113    #Channel of interest (not relevant for photonuclear processes)");
-    sl->SetParameter(Form("RND_SEED = %d  #Random number seed", sseed));
-    sl->SetParameter("OUTPUT_FORMAT =   2    #Form of the output");
-    sl->SetParameter("BREAKUP_MODE  =   5    #Controls the nuclear breakup");
-    sl->SetParameter("INTERFERENCE  =   0    #Interference (0 = off, 1 = on)");
-    sl->SetParameter("IF_STRENGTH   =   1.   #% of intefernce (0.0 - 0.1)");
-    sl->SetParameter("COHERENT      =   1    #Coherent=1,Incoherent=0");
-    sl->SetParameter("INCO_FACTOR   =   1.   #percentage of incoherence");
-    sl->SetParameter("BFORD         =   9.5  #");
-    sl->SetParameter("INT_PT_MAX    =   0.24 #Maximum pt considered, when interference is turned on");
-    sl->SetParameter("INT_PT_N_BINS = 120    #Number of pt bins when interference is turned on");
-
-    sl->Init();
-    sl->GetTStarLight()->PrintInputs(std::cout);
-    
-    AliGenerator* gener = sl;
-    gener->SetOrigin(0, 0, 0);    // vertex position
-    gener->SetSigma(0, 0, 5.3);   // Sigma in (X,Y,Z) (cm) on IP position
-    gener->SetCutVertexZ(1.);     // Truncate at 1 sigma
-    gener->SetVertexSmear(kPerEvent); 
-    gener->SetTrackingFlag(1);
-    gener->Init();
-    
-    if (smag == AliMagF::k2kG) {
-       comment = comment.Append(" | L3 field 0.2 T");
-    } else if (smag == AliMagF::k5kG) {
-       comment = comment.Append(" | L3 field 0.5 T");
-    }
-    
-    
-    printf("\n \n Comment: %s \n \n", comment.Data());
-    
-    
-// Field
-    TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", -1., -1., smag));
-
-    rl->CdGAFile();
-//
-    Int_t   iABSO   = 1;
-    Int_t   iDIPO   = 1;
-    Int_t   iFMD    = 1;
-    Int_t   iFRAME  = 1;
-    Int_t   iHALL   = 1;
-    Int_t   iITS    = 1;
-    Int_t   iMAG    = 1;
-    Int_t   iMUON   = 1;
-    Int_t   iPHOS   = 1;
-    Int_t   iPIPE   = 1;
-    Int_t   iPMD    = 1;
-    Int_t   iHMPID  = 1;
-    Int_t   iSHIL   = 1;
-    Int_t   iT0     = 1;
-    Int_t   iTOF    = 1;
-    Int_t   iTPC    = 1;
-    Int_t   iTRD    = 1;
-    Int_t   iZDC    = 1;
-    Int_t   iEMCAL  = 1;
-    Int_t   iVZERO  = 1;
-    Int_t   iACORDE = 1;
-
-    //=================== Alice BODY parameters =============================
-    AliBODY *BODY = new AliBODY("BODY", "Alice envelop");
-
-
-    if (iMAG)
-    {
-        //=================== MAG parameters ============================
-        // --- Start with Magnet since detector layouts may be depending ---
-        // --- on the selected Magnet dimensions ---
-        AliMAG *MAG = new AliMAG("MAG", "Magnet");
-    }
-
-
-    if (iABSO)
-    {
-        //=================== ABSO parameters ============================
-        AliABSO *ABSO = new AliABSOv3("ABSO", "Muon Absorber");
-    }
-
-    if (iDIPO)
-    {
-        //=================== DIPO parameters ============================
-
-        AliDIPO *DIPO = new AliDIPOv3("DIPO", "Dipole version 3");
-    }
-
-    if (iHALL)
-    {
-        //=================== HALL parameters ============================
-
-        AliHALL *HALL = new AliHALLv3("HALL", "Alice Hall");
-    }
-
-
-    if (iFRAME)
-    {
-        //=================== FRAME parameters ============================
-
-        AliFRAMEv2 *FRAME = new AliFRAMEv2("FRAME", "Space Frame");
-       FRAME->SetHoles(1);
-    }
-
-    if (iSHIL)
-    {
-        //=================== SHIL parameters ============================
-
-        AliSHIL *SHIL = new AliSHILv3("SHIL", "Shielding Version 3");
-    }
-
-
-    if (iPIPE)
-    {
-        //=================== PIPE parameters ============================
-
-        AliPIPE *PIPE = new AliPIPEv3("PIPE", "Beam Pipe");
-    }
-    if (iITS)
-    {
-        //=================== ITS parameters ============================
-
-       AliITS *ITS  = new AliITSv11("ITS","ITS v11");
-    }
-
-    if (iTPC)
-    {
-      //============================ TPC parameters =====================
-        AliTPC *TPC = new AliTPCv2("TPC", "Default");
-    }
-
-
-    if (iTOF) {
-        //=================== TOF parameters ============================
-       AliTOF *TOF = new AliTOFv6T0("TOF", "normal TOF");
-    }
-
-
-    if (iHMPID)
-    {
-        //=================== HMPID parameters ===========================
-        AliHMPID *HMPID = new AliHMPIDv3("HMPID", "normal HMPID");
-
-    }
-
-
-    if (iZDC)
-    {
-        //=================== ZDC parameters ============================
-
-        AliZDC *ZDC = new AliZDCv4("ZDC", "normal ZDC");
-    }
-
-    if (iTRD)
-    {
-        //=================== TRD parameters ============================
-
-        AliTRD *TRD = new AliTRDv1("TRD", "TRD slow simulator");
-    }
-
-    if (iFMD)
-    {
-        //=================== FMD parameters ============================
-       AliFMD *FMD = new AliFMDv1("FMD", "normal FMD");
-   }
-
-    if (iMUON)
-    {
-        //=================== MUON parameters ===========================
-        // New MUONv1 version (geometry defined via builders)
-        AliMUON *MUON = new AliMUONv1("MUON", "default");
-    }
-    //=================== PHOS parameters ===========================
-
-    if (iPHOS)
-    {
-        AliPHOS *PHOS = new AliPHOSv1("PHOS", "IHEP");
-    }
-
-
-    if (iPMD)
-    {
-        //=================== PMD parameters ============================
-        AliPMD *PMD = new AliPMDv1("PMD", "normal PMD");
-    }
-
-    if (iT0)
-    {
-        //=================== T0 parameters ============================
-        AliT0 *T0 = new AliT0v1("T0", "T0 Detector");
-    }
-
-    if (iEMCAL)
-    {
-        //=================== EMCAL parameters ============================
-        AliEMCAL *EMCAL = new AliEMCALv2("EMCAL", "EMCAL_COMPLETEV1");
-    }
-
-     if (iACORDE)
-    {
-        //=================== ACORDE parameters ============================
-        AliACORDE *ACORDE = new AliACORDEv1("ACORDE", "normal ACORDE");
-    }
-
-     if (iVZERO)
-    {
-        //=================== VZERO parameters ============================
-        AliVZERO *VZERO = new AliVZEROv7("VZERO", "normal VZERO");
-    }
-             
-}
-
-Float_t EtaToTheta(Float_t arg){
-  return (180./TMath::Pi())*2.*atan(exp(-arg));
-}
-
-
-void ProcessEnvironmentVars()
-{
-    // Run type
-    if (gSystem->Getenv("CONFIG_RUN_TYPE")) {
-      for (Int_t iRun = 0; iRun < kRunMax; iRun++) {
-       if (strcmp(gSystem->Getenv("CONFIG_RUN_TYPE"), pprRunName[iRun])==0) {
-         srun = (PprRun_t)iRun;
-         cout<<"Run type set to "<<pprRunName[iRun]<<endl;
-       }
-      }
-    }
-
-    // Random Number seed
-    if (gSystem->Getenv("CONFIG_SEED")) {
-      sseed = atoi(gSystem->Getenv("CONFIG_SEED"));
-    }
-}
diff --git a/STARLIGHT/test/galice.root b/STARLIGHT/test/galice.root
deleted file mode 100644 (file)
index eeeb037..0000000
Binary files a/STARLIGHT/test/galice.root and /dev/null differ
diff --git a/STARLIGHT/test/rec.C~ b/STARLIGHT/test/rec.C~
deleted file mode 100644 (file)
index bf8446d..0000000
+++ /dev/null
@@ -1,26 +0,0 @@
-void rec() {
-  AliReconstruction reco;
-
-  reco.SetWriteESDfriend();
-  reco.SetWriteAlignmentData();
-
-  reco.SetDefaultStorage("local://$ALICE_ROOT/OCDB");
-  reco.SetSpecificStorage("GRP/GRP/Data",
-                         Form("local://%s",gSystem->pwd()));
-  reco.SetRunPlaneEff(kTRUE);
-
-  reco.SetRunQA("ALL:ALL") ;
-  
-  reco.SetQARefDefaultStorage("local://$ALICE_ROOT/QAref") ;
-  
-  for (Int_t det = 0 ; det < AliQA::kNDET ; det++) {
-    reco.SetQACycles((AliQAv1::DETECTORINDEX_t)det, 999) ;
-    reco.SetQAWriteExpert((AliQAv1::DETECTORINDEX_t)det) ; 
-  }
-  
-  TStopwatch timer;
-  timer.Start();
-  reco.Run();
-  timer.Stop();
-  timer.Print();
-}
diff --git a/STARLIGHT/test/sim.C~ b/STARLIGHT/test/sim.C~
deleted file mode 100644 (file)
index baf8eeb..0000000
+++ /dev/null
@@ -1,36 +0,0 @@
-void sim(Int_t nev=20) {
-
-  gSystem->Load("liblhapdf");
-  gSystem->Load("libEGPythia6");
-  gSystem->Load("libpythia6");
-  gSystem->Load("libAliPythia6");
-  gSystem->Load("libhijing");
-  gSystem->Load("libTHijing");
-  gSystem->Load("libgeant321");
-
-  if (gSystem->Getenv("EVENT"))
-   nev = atoi(gSystem->Getenv("EVENT")) ;   
-  
-  AliSimulation simulator;
-  simulator.SetMakeSDigits("TRD TOF PHOS HMPID EMCAL MUON FMD ZDC PMD T0 VZERO");
-  simulator.SetMakeDigitsFromHits("ITS TPC");
-  simulator.SetWriteRawData("ALL","raw.root",kTRUE);
-
-  simulator.SetDefaultStorage("local://$ALICE_ROOT/OCDB");
-  simulator.SetSpecificStorage("GRP/GRP/Data",
-                              Form("local://%s",gSystem->pwd()));
-  
-  simulator.SetRunQA("ALL:ALL") ; 
-  
-  simulator.SetQARefDefaultStorage("local://$ALICE_ROOT/QAref") ;
-
-  for (Int_t det = 0 ; det < AliQA::kNDET ; det++) {
-    simulator.SetQACycles((AliQAv1::DETECTORINDEX_t)det, nev+1) ;
-  }
-  
-  TStopwatch timer;
-  timer.Start();
-  simulator.Run(nev);
-  timer.Stop();
-  timer.Print();
-}
diff --git a/STARLIGHT/test/syswatch.log b/STARLIGHT/test/syswatch.log
deleted file mode 100644 (file)
index 7adddd4..0000000
+++ /dev/null
@@ -1,6 +0,0 @@
-hname/C:sname/C:id0/I:id1/I:id2/I:id3/I:first/D:stampSec/D:mi.fMemUsed/D:mi.fSwapUsed/D:cI.fUser/D:cI.fSys/D:cI.fLoad1m/D:cI.fLoad5m/D:cI.fLoad15m/D:pI.fMemResident/D:pI.fMemVirtual/D:pI.fCpuUser/D:pI.fCpuSys/D:stampOldSec/D:miOld.fMemUsed/D:miOld.fSwapUsed/D:cIOld.fUser/D:cIOld.fSys/D:pIOld.fMemResident/D:pIOld.fMemVirtual/D:pIOld.fCpuUser/D:pIOld.fCpuSys/D:fileBytesRead/D:fileBytesWritten/D:fileCounter/D:fileReadCalls/D
-nz23-18.ifj.edu.pl     Start_Run       -1      -1      -1      -1      1385552498.02524        1385552498.02524        7859    192     100     0       1.02999997138977        1.03999996185303        0.920000016689301       113.0859375     224.42578125    0.422935009002686       0.084986999630928       1385552498.03631        0       0       0       0       0       0       0       0       0       0       0       0       
-nz23-18.ifj.edu.pl     ProcessEnvironmentVars  -1      -1      -1      -1      1385552498.02524        1385552498.03661        7861    192     50      50      1.02999997138977        1.03999996185303        0.920000016689301       113.1171875     224.42578125    0.423934996128082       0.084986999630928       1385552498.02524        7859    192     100     0       113.0859375     224.42578125    0.422935009002686       0.084986999630928       0       0       0       0       
-nz23-18.ifj.edu.pl     RunSimulation_Begin     -1      -1      -1      -1      1385552498.02524        1385552498.0568 7857    192     100     0       1.02999997138977        1.03999996185303        0.920000016689301       113.26953125    224.546875      0.428934007883072       0.0899859964847565      1385552498.03661        7861    192     50      50      113.1171875     224.42578125    0.423934996128082       0.084986999630928       0       0       0       0       
-nz23-18.ifj.edu.pl     RunSimulation_InitCDB   -1      -1      -1      -1      1385552498.02524        1385552499.99692        7864    192     100     0       1.02999997138977        1.03999996185303        0.920000016689301       119.953125      256.4453125     0.543917000293732       0.110982999205589       1385552498.0568 7857    192     100     0       113.26953125    224.546875      0.428934007883072       0.0899859964847565      0       0       0       0       
-nz23-18.ifj.edu.pl     RunSimulation_SetCDBLock        -1      -1      -1      -1      1385552498.02524        1385552501.99228        7863    192     100     0       1.02999997138977        1.03999996185303        0.920000016689301       120.41796875    256.7890625     0.550916016101837       0.110982999205589       1385552499.99692        7864    192     100     0       119.953125      256.4453125     0.543917000293732       0.110982999205589       0       0       0       0