PFireball.h

Go to the documentation of this file.

// Author: Romain Holzmann
// Written: 28.04.00
// Revised: 14.08.00 MK
// Revised: 05.09.05 RH
// Revised: 28.03.06 RH
// Revised: 02.02.07 RH
// Revised: 21.02.07 RH
// PFireball Class Header
 
#ifndef _PFireball_H_
#define _PFireball_H_
 
#include "PF2.h"
#include "TH2F.h"
#include "PParticle.h"
#include "PChannel.h"
#include "PChannelModel.h"
#include "PFunction.h"
#include "PAdaptiveMeshN.h"
 
class PFireball: public PParticle {
 
 protected:
 
    int prodId;    // product id
    float T1;      // temperature of component 1 (in GeV)
    float T2;      // temperature of component 2 (in GeV)
    float frac;    // fraction of component 1  (0 < frac < 1)
    float blast;   // radial blast velocity  (0 < blast < 1)
    float power1;  // power of p**n term in p**n exp(-E/T) sampling at 0 deg
    float power2;  // power of p**n term at 90 deg
    float A2;      // polar anisotropy (dN/dOmega = 1 + A2*cost**2 + A4*cost**4)
    float A4;      //
    float v1;      // directed flow parameter at mid-rapidity
    float v2;      // elliptic flow (squeeze out) at mid-rapidity
    float Ap;      // projectile mass
    float At;      // target mass
    float prob;    // particle production probability per participant nucleon
    float bmin;    // min impact parameter
    float bmax;    // max impact parameter
    float meanN;   // mean multiplicity
    int   nProd;   // number of product particles
    bool  flag;    
    float mt_fac;  
    double sig;    
    TF1  *rapidity_function; 
 
    int spect;     
    //  if set to 2 projectile-like spectator
 
    PF2 *fE;       
    TF1 *fE_1d;    
    PF2 *fE1;      
    TF1 *fE1_1d;   
    PF2 *fE2;      
    TF1 *fE2_1d;   
    PF2 *fE3;      
    TF1 *fA;       
    TF1 *fMt;      
    TH2F *fHisto;  
    bool trueThermal; 
    Int_t npx, npy; 
 
    PParticle **part; 
    bool quasistable; 
 
    PChannelModel *model;  
    int didx_old;          
    int sample_option;     
 
    int mesh_option;     
    PFunction *pfE;      
    PAdaptiveMeshN *afE; 
 
    Float_t y0;          // mid-rapidity of source
 
 public:
    PFireball(const char *particle, float AGeV, float t1, float t2=0.0, float f=1.0,
          float b=0.0, float a2=0.0, float a4=0.0, float w1=0.0, float w2=0.0, int sp=0);
 
    void setTemperature(float t1, float t2, float f, int id=0) {
    T1 = t1;
    T2 = t2;
    frac = f;
    updateFunctions(id);
    if (mt_fac > 0.)
        mt_fac = mtIntegral(makeStaticData()->GetParticleMass(prodId),T1);
    }
 
    //Rapidity sampling:
    void setSigma(double s) {
    SetSigma(s);
    }
    void SetSigma(double s) {
    sig = s;
    }
    void SetRapidityDistribution(TF1 *f) {
    rapidity_function = f;
    };
 
    void setBlast(float b, int id=0) {
    blast = b;
    updateFunctions(id);
    }
    void setAnisotropy(float a2, float a4, int id=0) {
    A2 = a2;
    A4 = a4;
    updateFunctions(id);
    }
    void setFlow(float a, float b) {
    v1 = a;
    v2 = b;
    }
 
    Double_t AvApart(Double_t ap, Double_t at, Double_t bn, Double_t bx);
    Double_t AvApart(Double_t ap, Double_t at) {
    return (ap*pow(at,0.667) + at*pow(ap,0.667))
        /pow((pow(ap,0.333)+pow(at,0.333)),2);  // average Apart
    }
    void setRandomB(float ap, float at, float p=0., float bn=0., float bx=20.) {
    Ap = ap;
    At = at;
    prob = p;
    flag = 1;
    bmin = bn;
    bmax = bx;
    if (prob <= 0.)
        prob = meanN/AvApart(Ap, At);
    meanN = prob*AvApart(Ap, At, bmin, bmax);
    if (bmax > 1.14*(pow((double)Ap,(double)0.333) + pow((double)At,(double)0.333))+2.)
        bmax = 1.14*(pow((double)Ap,(double)0.333) + pow((double)At,(double)0.333)+2.);
    if (bmin > bmax) {
        float temp = bmin;
        bmin = bmax;
        bmax = temp;
    }
    }
    void setMeanN(float mN) {
    if (mN >= 0.) meanN = mN;
    }
    void setSpectator(int sp){
    spect=sp;
    }
 
    bool IsRandomB() {return flag;}
    bool IsRandomN() {return meanN>0. ? 1 : 0;}
 
    int getParticleId() {return prodId;}
    float getT1()       {return T1;}
    float getT2()       {return T1;}
    float getFrac()     {return frac;}
    float getBlast()    {return blast;}
    float getA2()       {return A2;}
    float getA4()       {return A4;}
    float getV1()       {return v1;}
    float getV2()       {return v2;}
    int  getSpectator() {return spect;}
    PF2 *getFuncE()     {return fE;}
    PF2 *getFuncE1()    {return fE1;}
    PF2 *getFuncE2()    {return fE2;}
    PF2 *getFuncE3()    {return fE3;}
    TF1 *getFuncA()     {return fA;}
    TF1 *getFuncMt()    {return fMt;}
    void setHisto(TH2F *pH) {
    fHisto = pH;
    A2 = A4 = 0.;
    }
    TH2F *getHisto() {return fHisto;}
 
    void SetNpx(Int_t my_npx);
    void SetNpy(Int_t my_npy);
    void SetEpsilon(Double_t e);
 
    void RotateSamplingPlane(void) {
    sample_option = 1;
    updateFunctions();
    }
    void UseMesh(void) {
    mesh_option = 1;
    }
 
 
    void sampleECM(Double_t &E, Double_t &M, int didx);
 
    void sampleECM1(Double_t &E, Double_t &M, int didx) {
    if (fE1) {
        if (didx != didx_old) {
        didx_old = didx; //prevent random initialization each time
        if (didx < 0)
            fE1->SetParameter(5, -1.1);
        else
            fE1->SetParameter(5, (double) didx);
        }
        fE1->GetRandom2(E, M);
    } else if (fE1_1d) {
        E = fE1_1d->GetRandom();
    } else {
        Error("sampleECM1", "fE1 not found");
    }
    }
 
    void sampleECM2(Double_t &E, Double_t &M, int didx) {
    if (fE2) {
        if (didx != didx_old) {
        didx_old = didx; //prevent random initialization each time
        if (didx<0)
            fE2->SetParameter(5, -1.1);
        else
            fE2->SetParameter(5, (double) didx);
        }
        fE2->GetRandom2(E, M);
    } else if (fE2_1d) {
        E = fE2_1d->GetRandom();
    } else {
        Error("sampleECM2", "fE2 not found");
    }
    }
 
    void sampleECM3(Double_t &E, Double_t &Theta) {
    fE3->GetRandom2(E, Theta);
    }
 
    double sampleThetaCM() {return fA->GetRandom();}
    double sampleMt()      {return fMt->GetRandom();}
    void samplePartCM(double &px, double &py, double &pz, double &E, int didx);
    float sampleB() const ;
    int sampleNProd();
    int sampleNProd(float b);
    int getLastNProd() {return nProd;}
    int IsFireball()   {return 1;}
 
    PChannel *makeChannel(Int_t nMax, Float_t nAverage=0.);
 
    void setTrueThermal(Bool_t flag=kTRUE) {
    trueThermal = flag;
    }
    void setMtScaling() {
    mt_fac = mtIntegral(makeStaticData()->GetParticleMass(prodId),T1);
    }
    float mtScale(double m) {
    return mt_fac>0. ? mtIntegral(m,T1)/mt_fac : 1.;
    }
    virtual void Print(const Option_t *delme="") const;
    void printAverages() const ;
    virtual ~PFireball() {
    delete fE;
    delete fE_1d;
    delete fE1;
    delete fE1_1d;
    delete fE2;
    delete fE2_1d;
    delete fA;
    delete fMt;
    }
    float mtIntegral(double mass, float temperature);
 
    void SetW(double w=1.) {
    if (*(makeStaticData()->GetBatchValue("_system_weight_version"))) {
        Info("SetW", "Use old weighting method (pure chain based)");
    }
    *(makeStaticData()->GetBatchValue("_system_weight_version")) = 0.;
    PParticle::SetW(w);
    };
 
 protected:
    void updateFunctions(int id = 0);
    void setToMidrapidity(float AGeV);
    float Npar(float ap, float at, float b);
    float NparSmeared(float ap, float at, float b) const ;
 
    ClassDef(PFireball, 1) // Pluto Fireball Class
 
};
#endif // _PFireball_H_