diff --git a/CMakeLists.txt b/CMakeLists.txt index 8d03796..fe0edf4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -10,6 +10,8 @@ add_subdirectory(GeneratorCosmics) add_subdirectory(GeneratorParam) +add_subdirectory(GeneratorSlowNucleons) + if (DEFINED ENV{HIJING_ROOT}) add_subdirectory(THijing) endif (DEFINED ENV{HIJING_ROOT}) diff --git a/GeneratorSlowNucleons/CMakeLists.txt b/GeneratorSlowNucleons/CMakeLists.txt new file mode 100644 index 0000000..3a8e3c2 --- /dev/null +++ b/GeneratorSlowNucleons/CMakeLists.txt @@ -0,0 +1,45 @@ +cmake_minimum_required(VERSION 3.0 FATAL_ERROR) +project(GeneratorSlowNucleons) + +# You need to tell CMake where to find the ROOT installation. This can be done in a number of ways: +# - ROOT built with classic configure/make use the provided $ROOTSYS/etc/cmake/FindROOT.cmake +# - ROOT built with CMake. Add in CMAKE_PREFIX_PATH the installation prefix for ROOT +list(APPEND CMAKE_PREFIX_PATH $ENV{ROOTSYS}) + +#---Locate the ROOT package and defines a number of variables (e.g. ROOT_INCLUDE_DIRS) +find_package(ROOT REQUIRED COMPONENTS EG) + +#---Define useful ROOT functions and macros (e.g. ROOT_GENERATE_DICTIONARY) +include(${ROOT_USE_FILE}) + +include_directories(${CMAKE_CURRENT_SOURCE_DIR}/.) + +set(HEADERS GeneratorSlowNucleons.h SlowNucleonModel.h SlowNucleonModelExp.h) + +ROOT_GENERATE_DICTIONARY(G__GeneratorSlowNucleons ${HEADERS} LINKDEF GeneratorSlowNucleonsLinkDef.h) + +#---Create a shared library with geneated dictionary +add_library(GeneratorSlowNucleons SHARED GeneratorSlowNucleons.cxx SlowNucleonModel.cxx SlowNucleonModelExp.cxx G__GeneratorSlowNucleons.cxx) +target_link_libraries(GeneratorSlowNucleons ${ROOT_LIBRARIES}) + + +set_target_properties(GeneratorSlowNucleons + PROPERTIES + PUBLIC_HEADER "${HEADERS}" ) + + +install(TARGETS GeneratorSlowNucleons + LIBRARY DESTINATION lib + PUBLIC_HEADER DESTINATION include) + +if (${ROOT_VERSION} VERSION_GREATER "6.0") + install( + FILES + ${CMAKE_CURRENT_BINARY_DIR}/libGeneratorSlowNucleons_rdict.pcm + ${CMAKE_CURRENT_BINARY_DIR}/libGeneratorSlowNucleons.rootmap + DESTINATION lib) +endif (${ROOT_VERSION} VERSION_GREATER "6.0") + +if(${CMAKE_SYSTEM} MATCHES Darwin) + set_target_properties(GeneratorSlowNucleons PROPERTIES LINK_FLAGS "-undefined dynamic_lookup") +endif(${CMAKE_SYSTEM} MATCHES Darwin) diff --git a/GeneratorSlowNucleons/GeneratorSlowNucleons.cxx b/GeneratorSlowNucleons/GeneratorSlowNucleons.cxx new file mode 100644 index 0000000..80d2dad --- /dev/null +++ b/GeneratorSlowNucleons/GeneratorSlowNucleons.cxx @@ -0,0 +1,428 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// +// Generator for slow nucleons in pA interactions. +// Source is modelled by a relativistic Maxwell distributions. +// In this case the number of slow nucleons is determined from the number of +// wounded nuclei using a realisation of AliSlowNucleonModel. Original code by +// Ferenc Sikler +// + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "GeneratorSlowNucleons.h" +#include "SlowNucleonModel.h" + +ClassImp(GeneratorSlowNucleons) + + GeneratorSlowNucleons::GeneratorSlowNucleons() + : TGenerator("GeneratorParam", "GeneratorParam"), fCMS(0.), fMomentum(0.), + fBeta(0.), fPmax(0.), fCharge(0), fProtonDirection(1.), fTemperatureG(0.), + fBetaSourceG(0.), fTemperatureB(0.), fBetaSourceB(0.), fNgp(0), fNgn(0), + fNbp(0), fNbn(0), fDebug(0), fDebugHist1(0), fDebugHist2(0), + fThetaDistribution(), fCosThetaGrayHist(), fCosTheta(), + fBeamCrossingAngle(0.), fBeamDivergence(0.), fBeamDivEvent(0.), + fSmearMode(2), fNcoll(0), fSlowNucleonModel(0) { + // Default constructor +} + +GeneratorSlowNucleons::GeneratorSlowNucleons(Int_t npart) + : TGenerator("GeneratorParam", "GeneratorParam"), fCMS(14000.), + fMomentum(0.), fBeta(0.), fPmax(10.), fCharge(1), fProtonDirection(1.), + fTemperatureG(0.05), fBetaSourceG(0.05), fTemperatureB(0.005), + fBetaSourceB(0.), fNgp(0), fNgn(0), fNbp(0), fNbn(0), fDebug(0), + fDebugHist1(0), fDebugHist2(0), fThetaDistribution(), fCosThetaGrayHist(), + fCosTheta(), fBeamCrossingAngle(0.), fBeamDivergence(0.), + fBeamDivEvent(0.), fSmearMode(2), fNcoll(0), + fSlowNucleonModel(new SlowNucleonModel()) + +{ + // Constructor + fName = "SlowNucleons"; + fTitle = "Generator for gray particles in pA collisions"; +} + +//____________________________________________________________ +GeneratorSlowNucleons::~GeneratorSlowNucleons() { + // Destructor + delete fSlowNucleonModel; +} + +void GeneratorSlowNucleons::SetProtonDirection(Float_t dir) { + // Set direction of the proton to change between pA (1) and Ap (-1) + fProtonDirection = dir / TMath::Abs(dir); +} + +void GeneratorSlowNucleons::Init() { + // + // Initialization + // + Double_t kMass = TDatabasePDG::Instance()->GetParticle(kProton)->Mass(); + fMomentum = fCMS / 2. * Float_t(fZTarget) / Float_t(fATarget); + fBeta = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum); + printf(" fMomentum %f fBeta %1.10f\n", fMomentum, fBeta); + if (fDebug) { + fDebugHist1 = + new TH2F("DebugHist1", "nu vs N_slow", 100, 0., 100., 20, 0., 20.); + fDebugHist2 = + new TH2F("DebugHist2", "b vs N_slow", 100, 0., 100., 15, 0., 15.); + fCosThetaGrayHist = + new TH1F("fCosThetaGrayHist", "Gray particles angle", 100, -1., 1.); + } + // + // non-uniform cos(theta) distribution + // + if (fThetaDistribution != 0) { + fCosTheta = + new TF1("fCosTheta", + "(2./3.14159265358979312)/(exp(2./3.14159265358979312)-exp(-2./" + "3.14159265358979312))*exp(2.*x/3.14159265358979312)", + -1., 1.); + } + + if (TMath::Abs(fBeamCrossingAngle) > 0.) + printf("\n GeneratorSlowNucleons: applying crossing angle %f mrad to slow " + "nucleons\n", + fBeamCrossingAngle * 1000.); +} + +void GeneratorSlowNucleons::FinishRun() { + // End of run action + // Show histogram for debugging if requested. + if (fDebug) { + TCanvas *c = new TCanvas("c", "Canvas 1", 400, 10, 600, 700); + c->Divide(2, 1); + c->cd(1); + fDebugHist1->Draw("colz"); + c->cd(2); + fDebugHist2->Draw(); + c->cd(3); + fCosThetaGrayHist->Draw(); + } +} + +void GeneratorSlowNucleons::GenerateEvent() { + // + // Generate one event + // + const Float_t mp = (TDatabasePDG::Instance()->GetParticle(kProton))->Mass(); + const Float_t mn = (TDatabasePDG::Instance()->GetParticle(kNeutron))->Mass(); + const Float_t k2PI = 2. * TMath::Pi(); + const Float_t kRaddeg = 180. / TMath::Pi(); + const Float_t kDegrad = TMath::Pi() / 180.; + fParticles->Clear(); + // + // Communication with Gray Particle Model + // + // if (fCollisionGeometry) { + // Float_t b = fCollisionGeometry->b; + + // (1) Sikler' model + if (fSmearMode == 0) + fSlowNucleonModel->GetNumberOfSlowNucleons(fNcoll, fNgp, fNgn, fNbp, fNbn); + // (2) Model inspired on exp. data at lower energy (Gallio-Oppedisano) + // --- smearing the Ncoll fron generator used as input + else if (fSmearMode == 1) + fSlowNucleonModel->GetNumberOfSlowNucleons2(fNcoll, fNgp, fNgn, fNbp, fNbn); + // --- smearing directly Nslow + else if (fSmearMode == 2) + fSlowNucleonModel->GetNumberOfSlowNucleons2s(fNcoll, fNgp, fNgn, fNbp, + fNbn); + if (fDebug) { + // printf("Collision Geometry %f %d %d %d %d\n", b, nn, nwn, nnw, nwnw); + printf("Slow nucleons: %d grayp %d grayn %d blackp %d blackn \n", fNgp, + fNgn, fNbp, fNbn); + fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fNcoll, 1.); + // fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.); + } + // + Float_t p[3] = {0., 0., 0.}, theta = 0; + Float_t polar[3] = {0., 0., 0.}; + Int_t nt, i, j; + Int_t kf; + Double_t energy; + + // Extracting 1 value per event for the divergence angle + Double_t rvec = gRandom->Gaus(0.0, 1.0); + fBeamDivEvent = fBeamDivergence * TMath::Abs(rvec); + if (TMath::Abs(fBeamDivEvent) > 0.) + printf("\n GeneratorSlowNucleons: applying beam divergence %f mrad to " + "slow nucleons\n", + fBeamDivEvent * 1000.); + // + // Gray protons + // + fCharge = 1; + kf = kProton; + for (i = 0; i < fNgp; i++) { + GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta); + if (fDebug) + fCosThetaGrayHist->Fill(TMath::Cos(theta)); + energy = TMath::Sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2] + mp * mp); + auto part = new TParticle(kf, 1, -1, -1, -1, -1, p[0], p[1], p[2], energy, + 0., 0., 0., 0.); + part->SetUniqueID(kGrayProcess); + fParticles->Add(part); + } + // + // Gray neutrons + // + fCharge = 0; + kf = kNeutron; + for (i = 0; i < fNgn; i++) { + GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta); + if (fDebug) + fCosThetaGrayHist->Fill(TMath::Cos(theta)); + energy = TMath::Sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2] + mn * mn); + auto part = new TParticle(kf, 1, -1, -1, -1, -1, p[0], p[1], p[2], energy, + 0., 0., 0., 0.); + part->SetUniqueID(kGrayProcess); + fParticles->Add(part); + } + // + // Black protons + // + fCharge = 1; + kf = kProton; + for (i = 0; i < fNbp; i++) { + GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta); + energy = TMath::Sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2] + mp * mp); + auto part = new TParticle(kf, 1, -1, -1, -1, -1, p[0], p[1], p[2], energy, + 0., 0., 0., 0.); + part->SetUniqueID(kBlackProcess); + fParticles->Add(part); + } + // + // Black neutrons + // + fCharge = 0; + kf = kNeutron; + for (i = 0; i < fNbn; i++) { + GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta); + energy = TMath::Sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2] + mn * mn); + auto part = new TParticle(kf, 1, -1, -1, -1, -1, p[0], p[1], p[2], energy, + 0., 0., 0., 0.); + part->SetUniqueID(kBlackProcess); + fParticles->Add(part); + } +} + +int GeneratorSlowNucleons::ImportParticles(TClonesArray *particles, + Option_t *option) { + if (particles == 0) + return 0; + TClonesArray &clonesParticles = *particles; + clonesParticles.Clear(); + Int_t numpart = fParticles->GetEntries(); + for (int i = 0; i < numpart; i++) { + TParticle *particle = (TParticle *)fParticles->At(i); + new (clonesParticles[i]) TParticle(*particle); + } + return numpart; +} + +void GeneratorSlowNucleons::GenerateSlow(Int_t charge, Double_t T, + Double_t beta, Float_t *q, + Float_t &theta) + +{ + /* + Emit a slow nucleon with "temperature" T [GeV], + from a source moving with velocity beta + Three-momentum [GeV/c] is given back in q[3] + */ + + // printf("Generating slow nuc. with: charge %d. temp. %1.4f, beta %f + // \n",charge,T,beta); + + Double_t m, pmax, p, f, phi; + TDatabasePDG *pdg = TDatabasePDG::Instance(); + const Double_t kMassProton = pdg->GetParticle(kProton)->Mass(); + const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass(); + + /* Select nucleon type */ + if (charge == 0) + m = kMassNeutron; + else + m = kMassProton; + + /* Momentum at maximum of Maxwell-distribution */ + + pmax = TMath::Sqrt(2. * T * (T + TMath::Sqrt(T * T + m * m))); + + /* Try until proper momentum */ + /* for lack of primitive function of the Maxwell-distribution */ + /* a brute force trial-accept loop, normalized at pmax */ + + do { + p = gRandom->Rndm() * fPmax; + f = Maxwell(m, p, T) / Maxwell(m, pmax, T); + } while (f < gRandom->Rndm()); + + /* Spherical symmetric emission for black particles (beta=0)*/ + if (beta == 0 || fThetaDistribution == 0) + theta = TMath::ACos(2. * gRandom->Rndm() - 1.); + /* cos theta distributed according to experimental results for gray particles + * (beta=0.05)*/ + else if (fThetaDistribution != 0) { + Double_t costheta = fCosTheta->GetRandom(); + theta = TMath::ACos(costheta); + } + // + phi = 2. * TMath::Pi() * gRandom->Rndm(); + + /* Determine momentum components in system of the moving source */ + q[0] = p * TMath::Sin(theta) * TMath::Cos(phi); + q[1] = p * TMath::Sin(theta) * TMath::Sin(phi); + q[2] = p * TMath::Cos(theta); + if (fDebug == 1) + printf("\n Momentum in RS of the moving source: p = (%f, %f, %f)\n", q[0], + q[1], q[2]); + + /* Transform to system of the target nucleus */ + /* beta is passed as negative, because the gray nucleons are slowed down */ + Lorentz(m, -beta, q); + if (fDebug == 1) + printf(" Momentum in RS of the target nucleus: p = (%f, %f, %f)\n", q[0], + q[1], q[2]); + + /* Transform to laboratory system */ + Lorentz(m, fBeta, q); + q[2] *= fProtonDirection; + if (fDebug == 1) + printf("\n Momentum after LHC boost: p = (%f, %f, %f)\n", q[0], q[1], q[2]); + + if (TMath::Abs(fBeamCrossingAngle) > 0.) + BeamCrossing(q); // applying crossing angle + if (TMath::Abs(fBeamDivergence) > 0.) + BeamDivergence(q); // applying divergence +} + +Double_t GeneratorSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T) { + /* Relativistic Maxwell-distribution */ + Double_t ekin; + ekin = TMath::Sqrt(p * p + m * m) - m; + return (p * p * exp(-ekin / T)); +} + +//_____________________________________________________________________________ +void GeneratorSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t *q) { + /* Lorentz transform in the direction of q[2] */ + + Double_t gamma = 1. / TMath::Sqrt((1. - beta) * (1. + beta)); + Double_t energy = + TMath::Sqrt(m * m + q[0] * q[0] + q[1] * q[1] + q[2] * q[2]); + q[2] = gamma * (q[2] + beta * energy); +} + +//_____________________________________________________________________________ +void GeneratorSlowNucleons::BeamCrossing(Float_t *pLab) { + // Applying beam crossing angle + // + pLab[1] = pLab[2] * TMath::Sin(fBeamCrossingAngle) + + pLab[1] * TMath::Cos(fBeamCrossingAngle); + pLab[2] = pLab[2] * TMath::Cos(fBeamCrossingAngle) - + pLab[1] * TMath::Sin(fBeamCrossingAngle); + if (fDebug == 1) { + printf(" Beam crossing angle = %f mrad ", fBeamCrossingAngle * 1000.); + printf(" p = (%f, %f, %f)\n", pLab[0], pLab[1], pLab[2]); + } +} +//_____________________________________________________________________________ +void GeneratorSlowNucleons::BeamDivergence(Float_t *pLab) { + // Applying beam divergence + // + const Float_t k2PI = 2. * TMath::Pi(); + const Float_t kRaddeg = 180. / TMath::Pi(); + const Float_t kDegrad = TMath::Pi() / 180.; + Double_t pmod = + TMath::Sqrt(pLab[0] * pLab[0] + pLab[1] * pLab[1] + pLab[2] * pLab[2]); + + Double_t tetdiv = fBeamDivEvent; + Double_t fidiv = (gRandom->Rndm()) * 2. * TMath::Pi(); + + Double_t tetpart = + TMath::ATan2(TMath::Sqrt(pLab[0] * pLab[0] + pLab[1] * pLab[1]), pLab[2]); + Double_t fipart = 0.; + if (TMath::Abs(pLab[1]) > 0. || TMath::Abs(pLab[0]) > 0.) + fipart = TMath::ATan2(pLab[1], pLab[0]); + if (fipart < 0.) { + fipart = fipart + 2. * TMath::Pi(); + } + tetdiv = tetdiv * kRaddeg; + fidiv = fidiv * kRaddeg; + tetpart = tetpart * kRaddeg; + fipart = fipart * kRaddeg; + + Double_t angleSum[2] = {0., 0.}; + AddAngle(tetpart, fipart, tetdiv, fidiv, angleSum); + + Double_t tetsum = angleSum[0]; + Double_t fisum = angleSum[1]; + // printf("tetpart %f fipart %f tetdiv %f fidiv %f angleSum %f + // %f\n",tetpart,fipart,tetdiv,fidiv,angleSum[0],angleSum[1]); + tetsum = tetsum * kDegrad; + fisum = fisum * kDegrad; + + pLab[0] = pmod * TMath::Sin(tetsum) * TMath::Cos(fisum); + pLab[1] = pmod * TMath::Sin(tetsum) * TMath::Sin(fisum); + pLab[2] = pmod * TMath::Cos(tetsum); + if (fDebug == 1) { + printf(" Beam divergence %f mrad ", fBeamDivEvent * 1000.); + printf(" p = (%f, %f, %f)\n", pLab[0], pLab[1], pLab[2]); + } +} + +//_____________________________________________________________________________ +void GeneratorSlowNucleons::AddAngle(Double_t theta1, Double_t phi1, + Double_t theta2, Double_t phi2, + Double_t *angleSum) { + // Calculating the sum of 2 angles + Double_t temp = -1.; + Double_t conv = 180. / TMath::ACos(temp); + + Double_t ct1 = TMath::Cos(theta1 / conv); + Double_t st1 = TMath::Sin(theta1 / conv); + Double_t cp1 = TMath::Cos(phi1 / conv); + Double_t sp1 = TMath::Sin(phi1 / conv); + Double_t ct2 = TMath::Cos(theta2 / conv); + Double_t st2 = TMath::Sin(theta2 / conv); + Double_t cp2 = TMath::Cos(phi2 / conv); + Double_t sp2 = TMath::Sin(phi2 / conv); + Double_t cx = ct1 * cp1 * st2 * cp2 + st1 * cp1 * ct2 - sp1 * st2 * sp2; + Double_t cy = ct1 * sp1 * st2 * cp2 + st1 * sp1 * ct2 + cp1 * st2 * sp2; + Double_t cz = ct1 * ct2 - st1 * st2 * cp2; + + Double_t rtetsum = TMath::ACos(cz); + Double_t tetsum = conv * rtetsum; + if (TMath::Abs(tetsum) < 1e-4 || tetsum == 180.) + return; + + temp = cx / TMath::Sin(rtetsum); + if (temp > 1.) + temp = 1.; + if (temp < -1.) + temp = -1.; + Double_t fisum = conv * TMath::ACos(temp); + if (cy < 0) { + fisum = 360. - fisum; + } + + angleSum[0] = tetsum; + angleSum[1] = fisum; +} diff --git a/GeneratorSlowNucleons/GeneratorSlowNucleons.h b/GeneratorSlowNucleons/GeneratorSlowNucleons.h new file mode 100644 index 0000000..6dbadd7 --- /dev/null +++ b/GeneratorSlowNucleons/GeneratorSlowNucleons.h @@ -0,0 +1,137 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef O2_GENERATORSLOWNUCLEONS_H +#define O2_GENERATORSLOWNUCLEONS_H + +// +// Generator for slow nucleons in pA interactions. +// Source is modelled by a relativistic Maxwell distributions. +// Original code by Ferenc Sikler +// This class: andreas.morsch@cern.ch +// +#include + +class SlowNucleonModel; +class TH2F; +class TH1F; +class TF1; + +class GeneratorSlowNucleons : public TGenerator { +public: + GeneratorSlowNucleons(); + GeneratorSlowNucleons(Int_t npart); + virtual ~GeneratorSlowNucleons(); + virtual void Init(); + virtual void GenerateEvent(); + virtual int ImportParticles(TClonesArray *particles, Option_t *option); + virtual void FinishRun(); + virtual void SetPmax(Float_t pmax = 10.) { fPmax = pmax; } + virtual void SetNominalCmsEnergy(Float_t energy = 14000.) { fCMS = energy; } + virtual void SetTarget(Int_t a = 208, Int_t z = 82) { + fATarget = a; + fZTarget = z; + } + // virtual void SetTarget(TString s, Int_t a, Int_t z) + // {AliGenerator::SetTarget(s, a, z);} + virtual void SetProtonDirection(Float_t dir = 1.); + virtual void SetCharge(Int_t c = 1) { fCharge = c; } + virtual void SetTemperature(Double_t t1 = 0.04, Double_t t2 = 0.004) { + fTemperatureG = t1; + fTemperatureB = t2; + } + virtual void SetBetaSource(Double_t b1 = 0.05, Double_t b2 = 0.) { + fBetaSourceG = b1; + fBetaSourceB = b2; + } + // + virtual void SetSlowNucleonModel(SlowNucleonModel *model) { + fSlowNucleonModel = model; + } + virtual void SetNcoll(Int_t ncoll) { fNcoll = ncoll; } + virtual void SetDebug(Int_t flag = 0) { fDebug = flag; } + virtual void SetNumbersOfSlowNucleons(Int_t ngp, Int_t ngn, Int_t nbp, + Int_t nbn) { + fNgp = ngp; + fNgn = ngn; + fNbp = nbp; + fNbn = nbn; + } + // + // Added by Chiara to take into account angular distribution 4 gray tracks + virtual void SetThetaDist(Int_t flag = 0) { fThetaDistribution = flag; } + // + virtual void SetBeamCrossingAngle(Float_t crossAngle) { + fBeamCrossingAngle = crossAngle; + } + virtual void SetBeamDivergence(Float_t divergence) { + fBeamDivergence = divergence; + } + // + virtual Int_t GetNGrayProtons() { return fNgp; } + virtual Int_t GetNGrayNeutrons() { return fNgn; } + virtual Int_t GetNBlackProtons() { return fNbp; } + virtual Int_t GetNBlackNeutrons() { return fNbn; } + // + virtual void SetModelSmear(Int_t imode) { fSmearMode = imode; } + +protected: + void GenerateSlow(Int_t charge, Double_t T, Double_t beta, Float_t *q, + Float_t &theta); + Double_t Maxwell(Double_t m, Double_t p, Double_t t); + void Lorentz(Double_t m, Double_t beta, Float_t *q); + void BeamDivergence(Float_t *pLab); + void BeamCrossing(Float_t *pLab); + void AddAngle(Double_t theta1, Double_t phi1, Double_t theta2, Double_t phi2, + Double_t *angle); + +protected: + Float_t fCMS; // Center of mass energy + Double_t fMomentum; // Target nucleus momentum + Double_t fBeta; // Target nucleus beta + Float_t fPmax; // Maximum slow nucleon momentum + Int_t fCharge; // Slow nucleon charge + Int_t fATarget; // Target nucleus A + Int_t fZTarget; // Target nucleus Z + Float_t fProtonDirection; // Direction of the proton + Float_t fTemperatureG; // Source Temperature for gray nucleons + Float_t fBetaSourceG; // Source beta for gray nucleons + Float_t fTemperatureB; // Source Temperature for black nucleons + Float_t fBetaSourceB; // Source beta for black nucleons + Int_t fNgp; // Number of gray protons + Int_t fNgn; // Number of gray neutrons + Int_t fNbp; // Number of black protons + Int_t fNbn; // Number of black neutrons + Int_t fDebug; // Debug flag + TH2F *fDebugHist1; // Histogram for debugging + TH2F *fDebugHist2; // Histogram for debugging + // Added by Chiara to take into account angular distribution 4 gray tracks + Int_t fThetaDistribution; // 0 -> flat dist., 1 -> fwd. peaked distribution + TH1F *fCosThetaGrayHist; // Histogram for debugging + TF1 *fCosTheta; // Function for non-uniform cos(theta) distribution + // + Float_t fBeamCrossingAngle; // beam crossing angle (in radians) + Float_t fBeamDivergence; // beam divergence (in radians) + Float_t fBeamDivEvent; // beam divergence (in radians) + // + Int_t fSmearMode; // 0=Skler (no smear), =1 smearing Ncoll, =2 smearing Nslow + // + Int_t fNcoll; // number of collisions provided by external generator + SlowNucleonModel *fSlowNucleonModel; // The slow nucleon model + + enum { kGrayProcess = 200, kBlackProcess = 300 }; + +private: + GeneratorSlowNucleons(const GeneratorSlowNucleons &sn); + GeneratorSlowNucleons &operator=(const GeneratorSlowNucleons &rhs); + + ClassDef(GeneratorSlowNucleons, 1) // Slow Nucleon Generator +}; +#endif diff --git a/GeneratorSlowNucleons/GeneratorSlowNucleonsLinkDef.h b/GeneratorSlowNucleons/GeneratorSlowNucleonsLinkDef.h new file mode 100644 index 0000000..3544ec7 --- /dev/null +++ b/GeneratorSlowNucleons/GeneratorSlowNucleonsLinkDef.h @@ -0,0 +1,12 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class GeneratorSlowNucleons+; +#pragma link C++ class SlowNucleonModel+; +#pragma link C++ class SlowNucleonModelExp+; + + +#endif diff --git a/GeneratorSlowNucleons/SlowNucleonModel.cxx b/GeneratorSlowNucleons/SlowNucleonModel.cxx new file mode 100644 index 0000000..6ad9c11 --- /dev/null +++ b/GeneratorSlowNucleons/SlowNucleonModel.cxx @@ -0,0 +1,11 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +#include "SlowNucleonModel.h" +ClassImp(SlowNucleonModel) diff --git a/GeneratorSlowNucleons/SlowNucleonModel.h b/GeneratorSlowNucleons/SlowNucleonModel.h new file mode 100644 index 0000000..60c395f --- /dev/null +++ b/GeneratorSlowNucleons/SlowNucleonModel.h @@ -0,0 +1,38 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef O2_SLOWNUCLEONMODEL +#define O2_SLOWNUCLEONMODEL + +#include "TObject.h" +class SlowNucleonModel : public TObject { +public: + SlowNucleonModel() { ; } + virtual ~SlowNucleonModel() { ; } + virtual void GetNumberOfSlowNucleons(Int_t /*ncoll*/, Int_t & /*ngp*/, + Int_t & /*ngn*/, Int_t & /*nbp*/, + Int_t & /*nbn*/) const { + ; + } + virtual void GetNumberOfSlowNucleons2(Int_t /*ncoll*/, Int_t & /*ngp*/, + Int_t & /*ngn*/, Int_t & /*nbp*/, + Int_t & /*nbn*/) const { + ; + } + virtual void GetNumberOfSlowNucleons2s(Int_t /*ncoll*/, Int_t & /*ngp*/, + Int_t & /*ngn*/, Int_t & /*nbp*/, + Int_t & /*nbn*/) const { + ; + } + +protected: + ClassDef(SlowNucleonModel, 1) // Gray Particle Model +}; +#endif diff --git a/GeneratorSlowNucleons/SlowNucleonModelExp.cxx b/GeneratorSlowNucleons/SlowNucleonModelExp.cxx new file mode 100644 index 0000000..796db90 --- /dev/null +++ b/GeneratorSlowNucleons/SlowNucleonModelExp.cxx @@ -0,0 +1,253 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +// +// Experimental data inspired Gray Particle Model for p-Pb collisions +// The number of gray nucleons is proportional to the number of collisions. +// The number of black nucleons is proportional to the number of collisions +// Fluctuations are calculated from a binomial distribution. +// Author: A.Morsch +// + +#include "SlowNucleonModelExp.h" +#include +#include + +ClassImp(SlowNucleonModelExp) + + SlowNucleonModelExp::SlowNucleonModelExp() + : fP(82), fN(126), fAlphaGray(2.3), fAlphaBlack(3.6), + fApplySaturation(kTRUE), fnGraySaturation(15), fnBlackSaturation(28), + fLCPparam(0.585), fSigmaSmear(0.25) { + // + // Default constructor + // + // + fSlownparam[0] = 60.; + fSlownparam[1] = 469.2; + fSlownparam[2] = 8.762; + /*printf("\n\n ******** Initializing slow nucleon model with parameters:\n"); + printf(" \t alpha_{gray} %1.2f alpha_{black} %1.2f\n",fAlphaGray, + fAlphaBlack); printf(" \t SATURATION %d w. %d (gray) %d (black) + \n\n",fApplySaturation,fnGraySaturation,fnBlackSaturation); printf(" \t LCP + parameter %f Slown parameters = {%f, %f, + %f}\n\n",fLCPparam,fSlownparam[0],fSlownparam[1],fSlownparam[2]); */ +} + +void SlowNucleonModelExp::GetNumberOfSlowNucleons(Int_t ncoll, Int_t &ngp, + Int_t &ngn, Int_t &nbp, + Int_t &nbn) const { + // + // Return the number of black and gray nucleons + // + // Number of collisions + + Float_t nu = (Float_t)(ncoll); + + // Mean number of gray nucleons + + Float_t nGray = fAlphaGray * nu; + Float_t nGrayNeutrons = nGray * fN / (fN + fP); + Float_t nGrayProtons = nGray - nGrayNeutrons; + + // Mean number of black nucleons + Float_t nBlack = 0.; + if (!fApplySaturation || (fApplySaturation && nGray < fnGraySaturation)) + nBlack = fAlphaBlack * nu; + else if (fApplySaturation && nGray >= fnGraySaturation) + nBlack = fnBlackSaturation; + Float_t nBlackNeutrons = nBlack * 0.84; + Float_t nBlackProtons = nBlack - nBlackNeutrons; + + // Actual number (including fluctuations) from binomial distribution + Double_t p; + + // gray neutrons + p = nGrayNeutrons / fN; + ngn = gRandom->Binomial((Int_t)fN, p); + + // gray protons + p = nGrayProtons / fP; + ngp = gRandom->Binomial((Int_t)fP, p); + + // black neutrons + p = nBlackNeutrons / fN; + nbn = gRandom->Binomial((Int_t)fN, p); + + // black protons + p = nBlackProtons / fP; + nbp = gRandom->Binomial((Int_t)fP, p); +} + +void SlowNucleonModelExp::GetNumberOfSlowNucleons2(Int_t ncoll, Int_t &ngp, + Int_t &ngn, Int_t &nbp, + Int_t &nbn) const { + // + // Return the number of black and gray nucleons + // + // Number of collisions + + // based on E910 model + // ================================================================ + + Float_t nu = (Float_t)(ncoll); + // + // nu = nu+1.*gRandom->Rndm(); + nu = gRandom->Gaus(nu, 0.5); + if (nu < 0.) + nu = 0.; + // + Float_t poverpd = 0.843; + Float_t zAu2zPb = 82. / 79.; + Float_t nGrayp = (-0.27 + 0.63 * nu - 0.0008 * nu * nu) * poverpd * zAu2zPb; + + // gray protons + Double_t p; + p = nGrayp / fP; + ngp = gRandom->Binomial((Int_t)fP, p); + // ngp = gRandom->Gaus(nGrayp, TMath::Sqrt(fP*p*(1-p))); + if (nGrayp < 0.) + ngp = 0; + + // Float_t blackovergray = 3./7.;// from spallation + Float_t blackovergray = 0.65; // from COSY + Float_t nBlackp = blackovergray * nGrayp; + + // black protons + p = nBlackp / fP; + nbp = gRandom->Binomial((Int_t)fP, p); + // nbp = gRandom->Gaus(nBlackp, TMath::Sqrt(fP*p*(1-p))); + if (nBlackp < 0.) + nbp = 0; + + if (nu < 3.) { + nGrayp = -0.836 + 0.9112 * nu - 0.05381 * nu * nu; + nBlackp = blackovergray * nGrayp; + } + + // printf(" \t Using LCP parameter %f Slown parameters = {%f, %f, + // %f}\n\n",fLCPparam,fSlownparam[0],fSlownparam[1],fSlownparam[2]); + Float_t nGrayNeutrons = 0.; + Float_t nBlackNeutrons = 0.; + Float_t cp = (nGrayp + nBlackp) / fLCPparam; + + if (cp > 0.) { + Float_t nSlow = fSlownparam[0] + fSlownparam[1] / (-fSlownparam[2] - cp); + Float_t paramRetta = + fSlownparam[0] + fSlownparam[1] / (-fSlownparam[2] - 3); + if (cp < 3.) + nSlow = 0. + (paramRetta - 0.) / (3. - 0.) * (cp - 0.); + + nGrayNeutrons = nSlow * 0.1; + nBlackNeutrons = nSlow - nGrayNeutrons; + } else { + // Sikler "pasturato" (qui non entra mai!!!!) + nGrayNeutrons = 0.47 * fAlphaGray * nu; + nBlackNeutrons = 0.88 * fAlphaBlack * nu; + // printf("nslowp=0 -> ncoll = %1.0f -> ngrayn = %1.0f nblackn = %1.0f \n", + // nu, nGrayNeutrons, nBlackNeutrons); + } + + // gray neutrons + p = nGrayNeutrons / fN; + // ngn = gRandom->Binomial((Int_t) fN, p); + ngn = gRandom->Gaus(nGrayNeutrons, TMath::Sqrt(fN * p * (1 - p))); + + // black neutrons + p = nBlackNeutrons / fN; + // nbn = gRandom->Binomial((Int_t) fN, p); + nbn = gRandom->Gaus(nBlackNeutrons, TMath::Sqrt(fN * p * (1 - p))); +} + +void SlowNucleonModelExp::GetNumberOfSlowNucleons2s(Int_t ncoll, Int_t &ngp, + Int_t &ngn, Int_t &nbp, + Int_t &nbn) const { + // + // Return the number of black and gray nucleons + // + // Number of collisions + + // based on E910 model + // ================================================================ + + Float_t nu = (Float_t)(ncoll); + // + Float_t poverpd = 0.843; + Float_t zAu2zPb = 82. / 79.; + Float_t grayp = (-0.27 + 0.63 * nu - 0.0008 * nu * nu) * poverpd * zAu2zPb; + Float_t nGrayp = gRandom->Gaus(grayp, fSigmaSmear); + if (nGrayp < 0.) + nGrayp = 0.; + + // gray protons + Double_t p = 0.; + p = nGrayp / fP; + ngp = gRandom->Binomial((Int_t)fP, p); + // ngp = gRandom->Gaus(nGrayp, TMath::Sqrt(fP*p*(1-p))); + if (nGrayp < 0.) + ngp = 0; + + // Float_t blackovergray = 3./7.;// from spallation + Float_t blackovergray = 0.65; // from COSY + // Float_t blackp = blackovergray*grayp; + // Float_t nBlackp = gRandom->Gaus(nblackp, fSigmaSmear); + Float_t nBlackp = blackovergray * nGrayp; + if (nBlackp < 0.) + nBlackp = 0.; + + // black protons + p = nBlackp / fP; + nbp = gRandom->Binomial((Int_t)fP, p); + // nbp = gRandom->Gaus(nBlackp, TMath::Sqrt(fP*p*(1-p))); + if (nBlackp < 0.) + nbp = 0; + + Float_t nGrayNeutrons = 0.; + Float_t nBlackNeutrons = 0.; + Float_t cp = (nGrayp + nBlackp) / fLCPparam; + + if (cp > 0.) { + Float_t nSlow = fSlownparam[0] + fSlownparam[1] / (-fSlownparam[2] - cp); + + nGrayNeutrons = nSlow * 0.1; + nBlackNeutrons = nSlow - nGrayNeutrons; + } else { + // Sikler "pasturato" (qui non entra mai!!!!) + nGrayNeutrons = 0.47 * fAlphaGray * nu; + nBlackNeutrons = 0.88 * fAlphaBlack * nu; + // printf("nslowp=0 -> ncoll = %1.0f -> ngrayn = %1.0f nblackn = %1.0f \n", + // nu, nGrayNeutrons, nBlackNeutrons); + } + // + if (nGrayNeutrons < 0.) + nGrayNeutrons = 0.; + if (nBlackNeutrons < 0.) + nBlackNeutrons = 0.; + + // gray neutrons + p = nGrayNeutrons / fN; + // ngn = gRandom->Binomial((Int_t) fN, p); + ngn = gRandom->Gaus(nGrayNeutrons, TMath::Sqrt(fN * p * (1 - p))); + if (nGrayNeutrons < 0.) + ngn = 0; + + // black neutrons + p = nBlackNeutrons / fN; + // nbn = gRandom->Binomial((Int_t) fN, p); + nbn = gRandom->Gaus(nBlackNeutrons, TMath::Sqrt(fN * p * (1 - p))); + if (nBlackNeutrons < 0.) + nbn = 0; +} + +void SlowNucleonModelExp::SetParameters(Float_t alpha1, Float_t alpha2) { + // Set the model parameters + fAlphaGray = alpha1; + fAlphaBlack = alpha2; +} diff --git a/GeneratorSlowNucleons/SlowNucleonModelExp.h b/GeneratorSlowNucleons/SlowNucleonModelExp.h new file mode 100644 index 0000000..ab87c21 --- /dev/null +++ b/GeneratorSlowNucleons/SlowNucleonModelExp.h @@ -0,0 +1,68 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#ifndef O2_SLOWNUCLEONMODELEXP +#define O2_SLOWNUCLEONMODELEXP +// +// Experimental data inspired Gray Particle Model for p-Pb collisions +// Fluctuations are calculated from a binomial distribution. +// Author: A.Morsch +// + +#include "SlowNucleonModel.h" + +class SlowNucleonModelExp : public SlowNucleonModel { +public: + SlowNucleonModelExp(); + virtual ~SlowNucleonModelExp() { ; } + virtual void GetNumberOfSlowNucleons(Int_t ncoll, Int_t &ngp, Int_t &ngn, + Int_t &nbp, Int_t &nbn) const; + virtual void GetNumberOfSlowNucleons2(Int_t ncoll, Int_t &ngp, Int_t &ngn, + Int_t &nbp, Int_t &nbn) const; + virtual void GetNumberOfSlowNucleons2s(Int_t ncoll, Int_t &ngp, Int_t &ngn, + Int_t &nbp, Int_t &nbn) const; + // 1st model + virtual void SetParameters(Float_t alpha1, Float_t alpha2); + virtual void SetSaturation(Bool_t saturation) { + fApplySaturation = saturation; + } + virtual void SetSaturationParams(Int_t ngray = 15, Int_t nblack = 28) { + fnGraySaturation = ngray; + fnBlackSaturation = nblack; + } + // 2nd model + virtual void SetLCPparam(Float_t al) { fLCPparam = al; } + virtual void SetNslowParams(Float_t a, Float_t b, Float_t c) { + fSlownparam[0] = a; + fSlownparam[1] = b; + fSlownparam[2] = c; + } + +protected: + Float_t fP; // Number of protons in the target + Float_t fN; // Number of neutrons in the target + Float_t fAlphaGray; // Proportionality between gray particles and number of + // collisions + Float_t fAlphaBlack; // Proportionality between black particles and number of + // collisions + Bool_t fApplySaturation; // If true apply satoration to N_black vs. N_gray + Int_t fnGraySaturation; // N_gray value for N_black saturation + Int_t fnBlackSaturation; // N_black saturation value + // + // Adding parameters for 2nd model that can be tuned during config + Float_t fLCPparam; // parameter to calculate LCP from + Float_t fSlownparam[3]; // parameters to calculate from LCP + // + // Adding parameter to smear the number of slow nucleons + Float_t fSigmaSmear; + + ClassDef(SlowNucleonModelExp, 4) // Gray Particle Model (Experiment inspired) +}; +#endif diff --git a/GeneratorSlowNucleons/test.C b/GeneratorSlowNucleons/test.C new file mode 100644 index 0000000..37fdcb9 --- /dev/null +++ b/GeneratorSlowNucleons/test.C @@ -0,0 +1,20 @@ +void test() +{ + GeneratorSlowNucleons* slow = new GeneratorSlowNucleons(); + slow->SetSlowNucleonModel(new SlowNucleonModelExp()); + slow->SetNcoll(7); + slow->SetNominalCmsEnergy(13000); + slow->SetTarget(); + slow->Init(); + slow->GenerateEvent(); + + auto particles = new TClonesArray("TParticle", 10); + auto n = slow->ImportParticles(particles, "all"); + printf("Number of generated particles %5d \n", n); + + for (Int_t i = 0; i < n; i++) { + auto p = *(TParticle*) (particles->At(i)); + printf("%5d %5d %5d %5d %5d\n", i, p.GetPdgCode(), p.GetFirstMother(), p.GetFirstDaughter(), p.GetLastDaughter()); + + } +}