X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EVGEN%2FAliGenGeVSim.cxx;h=af1a8bf81d5a924a3eb0b3cea6fa01bbccda4484;hb=f564dd5a1c8891dcbf1537a0064343a362fc624a;hp=d39e278efc1a47d8d3ba4e6cab9b73eb79d0ca45;hpb=93ecfc91f2812d9b7d04e9cdbd9a3b5894e3f82c;p=u%2Fmrichter%2FAliRoot.git diff --git a/EVGEN/AliGenGeVSim.cxx b/EVGEN/AliGenGeVSim.cxx index d39e278efc1..af1a8bf81d5 100644 --- a/EVGEN/AliGenGeVSim.cxx +++ b/EVGEN/AliGenGeVSim.cxx @@ -1,5 +1,20 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +/* $Id$ */ -//////////////////////////////////////////////////////////////////////////////// // // AliGenGeVSim is a class implementing GeVSim event generator. // @@ -49,46 +64,66 @@ //////////////////////////////////////////////////////////////////////////////// -#include "Riostream.h" - -#include "TROOT.h" -#include "TCanvas.h" -#include "TParticle.h" -#include "TObjArray.h" -#include "TF1.h" -#include "TF2.h" -#include "TH1.h" -#include "TH2.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include -#include "AliRun.h" -#include "AliPDG.h" -#include "AliGenerator.h" -#include "AliGenGeVSim.h" #include "AliGeVSimParticle.h" +#include "AliGenGeVSim.h" #include "AliGenGeVSimEventHeader.h" +#include "AliGenerator.h" +#include "AliRun.h" -ClassImp(AliGenGeVSim); +ClassImp(AliGenGeVSim) ////////////////////////////////////////////////////////////////////////////////// -AliGenGeVSim::AliGenGeVSim() : AliGenerator(-1) { +AliGenGeVSim::AliGenGeVSim() : + AliGenerator(-1), + fModel(0), + fPsi(0), + fIsMultTotal(kTRUE), + fPtFormula(0), + fYFormula(0), + fPhiFormula(0), + fCurrentForm(0), + fPtYHist(0), + fPartTypes(0) +{ // // Default constructor // - fPsi = 0; - fIsMultTotal = kTRUE; - - //PH InitFormula(); for (Int_t i=0; i<4; i++) fPtYFormula[i] = 0; + for (Int_t i=0; i<2; i++) + fHist[i] = 0; } ////////////////////////////////////////////////////////////////////////////////// -AliGenGeVSim::AliGenGeVSim(Float_t psi, Bool_t isMultTotal) : AliGenerator(-1) { +AliGenGeVSim::AliGenGeVSim(Float_t psi, Bool_t isMultTotal) + : AliGenerator(-1), + fModel(0), + fPsi(psi), + fIsMultTotal(isMultTotal), + fPtFormula(0), + fYFormula(0), + fPhiFormula(0), + fCurrentForm(0), + fPtYHist(0), + fPartTypes(0) + { // // Standard Constructor. // @@ -109,12 +144,12 @@ AliGenGeVSim::AliGenGeVSim(Float_t psi, Bool_t isMultTotal) : AliGenerator(-1) { // checking consistancy if (psi < 0 || psi > 360 ) - Error ("AliGenGeVSim", "Reaction plane angle ( %d )out of range [0..360]", psi); + Error ("AliGenGeVSim", "Reaction plane angle ( %13.3f )out of range [0..360]", psi); fPsi = psi * TMath::Pi() / 180. ; fIsMultTotal = isMultTotal; - // initialization + // Initialization fPartTypes = new TObjArray(); InitFormula(); @@ -135,7 +170,7 @@ AliGenGeVSim::~AliGenGeVSim() { ////////////////////////////////////////////////////////////////////////////////// -Bool_t AliGenGeVSim::CheckPtYPhi(Float_t pt, Float_t y, Float_t phi) { +Bool_t AliGenGeVSim::CheckPtYPhi(Float_t pt, Float_t y, Float_t phi) const { // // private function used by Generate() // @@ -177,24 +212,92 @@ Bool_t AliGenGeVSim::CheckAcceptance(Float_t p[3]) { ////////////////////////////////////////////////////////////////////////////////// +// Deconvoluted Pt Y formula + +static Double_t aPtForm(Double_t * x, Double_t * par) { + // ptForm: pt -> x[0] , mass -> [0] , temperature -> [1] + // Description as string: " x * exp( -sqrt([0]*[0] + x*x) / [1] )" + + return x[0] * TMath::Exp( -sqrt(par[0]*par[0] + x[0]*x[0]) / par[1]); + } + +static Double_t aYForm(Double_t * x, Double_t * par) { + // y Form: y -> x[0] , sigmaY -> [0] + // Description as string: " exp ( - x*x / (2 * [0]*[0] ) )" + + return TMath::Exp ( - x[0]*x[0] / (2 * par[0]*par[0] ) ); + } + +// Models 1-3 +// Description as strings: + +// const char *kFormE = " ( sqrt([0]*[0] + x*x) * cosh(y) ) "; +// const char *kFormG = " ( 1 / sqrt( 1 - [2]*[2] ) ) "; +// const char *kFormYp = "( [2]*sqrt(([0]*[0]+x*x)*cosh(y)*cosh(y)-[0]*[0])/([1]*sqrt(1-[2]*[2]))) "; + +// const char* kFormula[3] = { +// " x * %s * exp( -%s / [1]) ", +// " (x * %s) / ( exp( %s / [1]) - 1 ) ", +// " x*%s*exp(-%s*%s/[1])*((sinh(%s)/%s)+([1]/(%s*%s))*(sinh(%s)/%s-cosh(%s)))" +// }; +// printf(kFormula[0], kFormE, kFormE); +// printf(kFormula[1], kFormE, kFormE); +// printf(kFormula[2], kFormE, kFormG, kFormE, kFormYp, kFormYp, kFormG, kFormE, kFormYp, kFormYp, kFormYp); + + +static Double_t aPtYFormula0(Double_t *x, Double_t * par) { + // pt -> x , Y -> y + // mass -> [0] , temperature -> [1] , expansion velocity -> [2] + + Double_t aFormE = TMath::Sqrt(par[0]*par[0] + x[0]*x[0]) * TMath::CosH(x[1]); + return x[0] * aFormE * TMath::Exp(-aFormE/par[1]); +} + +static Double_t aPtYFormula1(Double_t *x, Double_t * par) { + // pt -> x , Y -> y + // mass -> [0] , temperature -> [1] , expansion velocity -> [2] + + Double_t aFormE = TMath::Sqrt(par[0]*par[0] + x[0]*x[0]) * TMath::CosH(x[1]); + return x[0] * aFormE / ( TMath::Exp( aFormE / par[1]) - 1 ); +} + +static Double_t aPtYFormula2(Double_t *x, Double_t * par) { + // pt -> x , Y -> y + // mass -> [0] , temperature -> [1] , expansion velocity -> [2] + + Double_t aFormE = TMath::Sqrt(par[0]*par[0] + x[0]*x[0]) * TMath::CosH(x[1]); + Double_t aFormG = 1 / TMath::Sqrt((1.-par[2])*(1.+par[2])); + Double_t aFormYp = par[2]*TMath::Sqrt( (par[0]*par[0] + x[0]*x[0]) + * (TMath::CosH(x[1])-par[0])*(TMath::CosH(x[1])+par[0])) + /( par[1]*TMath::Sqrt((1.-par[2])*(1.+par[2]))); + + return x[0] * aFormE * TMath::Exp( - aFormG * aFormE / par[1]) + *( TMath::SinH(aFormYp)/aFormYp + + par[1]/(aFormG*aFormE) + * ( TMath::SinH(aFormYp)/aFormYp-TMath::CosH(aFormYp) ) ); +} + +// Phi Flow Formula + +static Double_t aPhiForm(Double_t * x, Double_t * par) { + // phi -> x + // Psi -> [0] , Direct Flow -> [1] , Elliptical Flow -> [2] + // Description as string: " 1 + 2*[1]*cos(x-[0]) + 2*[2]*cos(2*(x-[0])) " + + return 1 + 2*par[1]*TMath::Cos(x[0]-par[0]) + + 2*par[2]*TMath::Cos(2*(x[0]-par[0])); +} + void AliGenGeVSim::InitFormula() { // // private function // // Initalizes formulas used in GeVSim. - // Manages strings and creates TFormula objects from strings - // // Deconvoluted Pt Y formula - // ptForm: pt -> x , mass -> [0] , temperature -> [1] - // y Form: y -> x , sigmaY -> [0] - - const char* ptForm = " x * exp( -sqrt([0]*[0] + x*x) / [1] )"; - const char* yForm = " exp ( - x*x / (2 * [0]*[0] ) )"; - - fPtFormula = new TF1("gevsimPt", ptForm, 0, 3); - fYFormula = new TF1("gevsimRapidity", yForm, -3, 3); + fPtFormula = new TF1("gevsimPt", &aPtForm, 0, 3, 2); + fYFormula = new TF1("gevsimRapidity", &aYForm, -3, 3,1); fPtFormula->SetParNames("mass", "temperature"); fPtFormula->SetParameters(1., 1.); @@ -209,38 +312,19 @@ void AliGenGeVSim::InitFormula() { // Models 1-3 - // pt -> x , Y -> y - // mass -> [0] , temperature -> [1] , expansion velocity -> [2] - - - const char *formE = " ( sqrt([0]*[0] + x*x) * cosh(y) ) "; - const char *formG = " ( 1 / sqrt( 1 - [2]*[2] ) ) "; - const char *formYp = "( [2]*sqrt(([0]*[0]+x*x)*cosh(y)*cosh(y)-[0]*[0])/([1]*sqrt(1-[2]*[2]))) "; - - const char* formula[3] = { - " x * %s * exp( -%s / [1]) ", - " (x * %s) / ( exp( %s / [1]) - 1 ) ", - " x*%s*exp(-%s*%s/[1])*((sinh(%s)/%s)+([1]/(%s*%s))*(sinh(%s)/%s-cosh(%s)))" - }; + fPtYFormula[0] = new TF2("gevsimPtY_2", &aPtYFormula0, 0, 3, -2, 2, 2); - const char* paramNames[3] = {"mass", "temperature", "expVel"}; + fPtYFormula[1] = new TF2("gevsimPtY_3", &aPtYFormula1, 0, 3, -2, 2, 2); - char buffer[1024]; - - sprintf(buffer, formula[0], formE, formE); - fPtYFormula[0] = new TF2("gevsimPtY_2", buffer, 0, 3, -2, 2); - - sprintf(buffer, formula[1], formE, formE); - fPtYFormula[1] = new TF2("gevsimPtY_3", buffer, 0, 3, -2, 2); - - sprintf(buffer, formula[2], formE, formG, formE, formYp, formYp, formG, formE, formYp, formYp, formYp); - fPtYFormula[2] = new TF2("gevsimPtY_4", buffer, 0, 3, -2, 2); + fPtYFormula[2] = new TF2("gevsimPtY_4", &aPtYFormula2, 0, 3, -2, 2, 3); fPtYFormula[3] = 0; // setting names & initialisation + const char* kParamNames[3] = {"mass", "temperature", "expVel"}; + Int_t i, j; for (i=0; i<3; i++) { @@ -250,18 +334,14 @@ void AliGenGeVSim::InitFormula() { for (j=0; j<3; j++) { if ( i != 2 && j == 2 ) continue; // ExpVel - fPtYFormula[i]->SetParName(j, paramNames[j]); + fPtYFormula[i]->SetParName(j, kParamNames[j]); fPtYFormula[i]->SetParameter(j, 0.5); } } // Phi Flow Formula - // phi -> x - // Psi -> [0] , Direct Flow -> [1] , Ellipticla Flow -> [2] - - const char* phiForm = " 1 + 2*[1]*cos(x-[0]) + 2*[2]*cos(2*(x-[0])) "; - fPhiFormula = new TF1("gevsimPhi", phiForm, 0, 2*TMath::Pi()); + fPhiFormula = new TF1("gevsimPhi", &aPhiForm, 0, 2*TMath::Pi(), 3); fPhiFormula->SetParNames("psi", "directed", "elliptic"); fPhiFormula->SetParameters(0., 0., 0.); @@ -405,11 +485,11 @@ Float_t AliGenGeVSim::GetdNdYToTotal() { // Float_t integ, mag; - const Double_t maxPt = 3.0, maxY = 2.; + const Double_t kMaxPt = 3.0, kMaxY = 2.; if (fModel == 1) { - integ = fYFormula->Integral(-maxY, maxY); + integ = fYFormula->Integral(-kMaxY, kMaxY); mag = fYFormula->Eval(0); return integ/mag; } @@ -418,8 +498,8 @@ Float_t AliGenGeVSim::GetdNdYToTotal() { if (fModel > 1 && fModel < 6) { - integ = ((TF2*)fCurrentForm)->Integral(0,maxPt, -maxY, maxY); - mag = ((TF2*)fCurrentForm)->Integral(0, maxPt, -0.1, 0.1) / 0.2; + integ = ((TF2*)fCurrentForm)->Integral(0,kMaxPt, -kMaxY, kMaxY); + mag = ((TF2*)fCurrentForm)->Integral(0, kMaxPt, -0.1, 0.1) / 0.2; return integ/mag; } @@ -538,7 +618,7 @@ void AliGenGeVSim::SetFormula(Int_t pdg) { fPtYHist = (TH2D*)gROOT->FindObject(buff); } - if (!fPtYHist) Error(where, msg[4], pdg); + if (!fPtYHist) Error(where, msg[3], pdg); } } @@ -685,7 +765,7 @@ void AliGenGeVSim::Generate() { AliGeVSimParticle *partType; Int_t nType, nParticle, nParam; - const Int_t nParams = 6; + const Int_t kNParams = 6; // reaction plane determination and model DetermineReactionPlane(); @@ -707,7 +787,7 @@ void AliGenGeVSim::Generate() { // Evaluation of parameters - loop over parameters - for (nParam = 0; nParam < nParams; nParam++) { + for (nParam = 0; nParam < kNParams; nParam++) { paramScaler = FindScaler(nParam, pdg); @@ -790,7 +870,7 @@ void AliGenGeVSim::Generate() { // putting particle on the stack - SetTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, fTrackIt); + PushTrack(fTrackIt, kParent, pdg, p, orgin, polar, time, kPPrimary, id, fTrackIt); if (isMultTotal) nParticle++; } }