StFastGlauberMcMaker.h

/******************************************************************************
 * $Id: StFastGlauberMcMaker.h,v 1.4 2013/02/05 22:49:50 hmasui Exp $
 * $Log: StFastGlauberMcMaker.h,v $
 * Revision 1.4  2013/02/05 22:49:50  hmasui
 * Added Pb and Cu nuclei
 *
 * Revision 1.3  2012/04/25 05:22:42  hmasui
 * Added deformation parameters, higher order participant eccentricity with r^2 weight
 *
 * Revision 1.2  2010/11/20 19:03:11  hmasui
 * Mode mode flag to StCentralityMaker. Use STAR logger rather than STD iostream
 *
******************************************************************************/
//====================================================================================================
//  StFastGlauberMcMaker class 
//    Fast MC Glauber simulation
//
//   In order to run the MC glauber simulation, you can do
//   root4star -b
//   [0] .L doFastGlauberMcMaker.C
//   [1] doFastGlauberMcMaker("glauber.root", 1000, "AuAu", 200, "default", kFALSE);
//   where 1st argument is output ROOT file name, 2nd is number of events, 
//   3rd is the system, 4th is the center of mass energy
//   and 5th is type. The available type can be checked by
//     StFastGlauberMcMaker::Print("type");
//
//   - The 5th argument is the flag of deformation. kFALSE (or false) gives 
//     spherical nuclei (default), kTRUE (or true) gives deformed nuclei.
//     Deformation parameters are defined inside the class up to 4th order deformation.
//
//   You can see the debugging messages by
//     StFastGlauberMcMaker::DebugOn()
//   in case you need further checks of outputs
//
//----------------------------------------------------------------------------------------------------
//  Hiroshi Masui (HMasui@lbl.gov)
//====================================================================================================

#ifndef __StFastGlauberMcMaker_h__
#define __StFastGlauberMcMaker_h__

class TF1 ;
class TF2 ;
class TF3 ;
class TGraph ;
class TH1 ;
class TTree ;
class Nucleon ;
class StCentralityMaker ;
class StGlauberTree ;
#include <vector>
#include "TString.h"

//____________________________________________________________________________________________________
// Class StFastGlauberMcMaker: Fast Glauber Monte Carlo
class StFastGlauberMcMaker {
  public:
    // Default constructor (spherical Au+Au), and output file "fastglaubermc.root"
    StFastGlauberMcMaker();

    // Current available system:
    //     AuAu
    //
    // Current available type
    //     default
    //     large              Large R(+2%), small d(-10%)
    //     small              Small R(-2%), large d(+10%)
    //     largeXsec          Large inelastic NN cross section (+1mb)
    //     smallXsec          Small inelastic NN cross section (-1mb)
    //     gauss              Use gaussian collision profile
    StFastGlauberMcMaker(
        const TString outputFileName, // Output fileName
        const TString system, // System name (e.x. AuAu)
        const Double_t energy, // energy (GeV)
        const TString type,  // type (see above)
        const Bool_t isDeformed = kFALSE // Deformation flag
        );

    StFastGlauberMcMaker(
        const TString outputFileName,
        const UInt_t massNumber,  // Mass number of nucleus
        const Double_t radius,    // Radius of nucleus
        const Double_t skinDepth, // Skin depth of nucleus
        const Double_t beta2,     // 2nd order deformation parameter
        const Double_t beta4,     // 4th order deformation parameter
        const Double_t inelasticCrossSection, // Inelastic NN cross section
        const Double_t energy  // sqrt(sNN)
        ); 

    StFastGlauberMcMaker(
        const TString outputFileName,
        const UInt_t massNumberA,  // Mass number of nucleus for nucleus A
        const Double_t radiusA,    // Radius of nucleus for nucleus A
        const Double_t skinDepthA, // Skin depth of nucleus for nucleus A
        const Double_t beta2A,     // 2nd order deformation parameter for nucleus A
        const Double_t beta4A,     // 4th order deformation parameter for nucleus A
        const UInt_t massNumberB,  // Mass number of nucleus for nucleus B
        const Double_t radiusB,    // Radius of nucleus for nucleus B
        const Double_t skinDepthB, // Skin depth of nucleus for nucleus B
        const Double_t beta2B,     // 2nd order deformation parameter for nucleus B
        const Double_t beta4B,     // 4th order deformation parameter for nucleus B
        const Double_t inelasticCrossSection, // Inelastic NN cross section
        const Double_t energy  // sqrt(sNN)
        ); 

    virtual ~StFastGlauberMcMaker(); 

    void SetRepulsionDistance(const Double_t repulsionDistance);

    Int_t Make() ; 
    Int_t Run(const UInt_t nevents) ; 
    Int_t Finish() ; 

    void DoHardCoreSmearing() ; 

    void DoGaussianSmearing() ; 

    void DoHardCoreCollision() ; 
    void DoGaussianCollision() ; 

    //  option=type   print all available types
    void Print(const TString option="") const ;

    void DebugOn() ; 
    UInt_t Version() const ; 

  private:
    // Functions
    Int_t Clear() ; 

    Int_t InitTree() ;

    Int_t Init(
        const UInt_t massNumberA,  // Mass number of nucleus for nucleus A
        const Double_t radiusA,    // Radius of nucleus for nucleus A
        const Double_t skinDepthA, // Skin depth of nucleus for nucleus A
        const Double_t beta2A,     // 2nd order deformation parameter for nucleus A
        const Double_t beta4A,     // 4th order deformation parameter for nucleus A
        const UInt_t massNumberB,  // Mass number of nucleus for nucleus B
        const Double_t radiusB,    // Radius of nucleus for nucleus B
        const Double_t skinDepthB, // Skin depth of nucleus for nucleus B
        const Double_t beta2B,     // 2nd order deformation parameter for nucleus B
        const Double_t beta4B,      // 4th order deformation parameter for nucleus B
        const TString type = "default"
        );

    Int_t InitAuAu(const TString type) ;
    Int_t InitSmSm(const TString type) ;
    Int_t InitUU(const TString type) ;
    Int_t InitPbPb(const TString type) ;
    Int_t InitCuCu(const TString type) ;
    Int_t InitZrZr(const TString type, int Case) ;
    Int_t InitRuRu(const TString type, int Case) ;

    Bool_t IsCollision(Nucleon* nucleon0, Nucleon* nucleon1) const ;

    const Char_t* GetName(const UInt_t massNumber) const ; 
    Double_t GetInelasticNNCrossSection(const Double_t energy, const TString type) const ;
    void GetRThetaPhi(const UInt_t inucleus,
        Double_t& r, Double_t& theta, Double_t& phi) const ; 

    void Smearing(Double_t& r, Double_t& theta, Double_t& phi) const ; 

    //  
    //  TGraph::GetY()[0]    <r_T^2 * cos(n*phi)>
    //  TGraph::GetY()[1]    <r_T^2 * sin(n*phi)>
    //  TGraph::GetY()[2]    n * Psi = tan^{-1}(q_{y;n}, q_{x;n})
    //  TGraph::GetY()[3]    eps_{part;n}
    TGraph* GetParticipantEccentricity(const Double_t order, const Double_t sumx, const Double_t sumy,
        const Double_t sumw, const UInt_t weightId) const ;

    //____________________________________________________________________________________________________
    // Data members

    // For collision profile
    enum {
      mkHardCoreProfile = 0,
      mkGaussianProfile = 1
    };

    // Version control
    static const UInt_t mVersion ;

    const Double_t mEnergy ; 
    const TString mOutputFileName ; 
    Double_t mInelasticNNCrossSection ; 

    UInt_t mDebug ; 

    // Switches
    Double_t mRepulsionDistance ; 

    // NOTE: Only one smearing method can be applied
    //       Another option will be switched off when one turn on one of them
    Bool_t mHardCoreSmearing    ; 
    Bool_t mGaussianSmearing    ; 
    UInt_t mCollisionProfile    ; 

    UInt_t mMode ; 
    Bool_t mIsDeformed[2] ; 

    StCentralityMaker* mCentralityMaker ; 
    std::vector<Nucleon*> mNucleons[2] ; 
    TF1* mfWoodsSaxon[2]   ; 
    TF2* mfWoodsSaxon2D[2] ; 

    UInt_t mNeventsThrow  ; 
    UInt_t mNeventsAccept ; 
    TTree* mHeader        ; 

    // Output data in tree
    StGlauberTree* mGlauberTree ; 

    // QA histograms
    TH1* mhWoodsSaxon[4] ; 

    // Random number generator
    TF1* mfB        ; 
    TF3* mfHardCore ; 
    TF3* mfGaussian ; 

    ClassDef(StFastGlauberMcMaker, 1)
};

#endif