Adding GeneratorSpectator: generator for spectator protons (#9)

* Adding a test macro for GeneratorSpectators

* Set pmax in def. constructor to avoid infinite loop in generation + small fixes

Author: coppedis

GeneratorSlowNucleons/GeneratorSlowNucleons.cxx

@@ -32,7 +32,7 @@ ClassImp(GeneratorSlowNucleons)
 
     GeneratorSlowNucleons::GeneratorSlowNucleons()
     : TGenerator("GeneratorParam", "GeneratorParam"), fCMS(0.), fMomentum(0.),
-      fBeta(0.), fPmax(0.), fCharge(0), fProtonDirection(1.), fTemperatureG(0.),
+      fBeta(0.), fPmax(10.), 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(),
@@ -247,10 +247,9 @@ void GeneratorSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
      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);
+  //printf("Generating slow nuc. with: charge %d. temp. %1.4f, beta %f \n",charge,T,beta);
 
-  Double_t m, pmax, p, f, phi;
+  Double_t m=0, pmax=0, p=0, f=0, phi=0;
   TDatabasePDG *pdg = TDatabasePDG::Instance();
   const Double_t kMassProton = pdg->GetParticle(kProton)->Mass();
   const Double_t kMassNeutron = pdg->GetParticle(kNeutron)->Mass();
@@ -315,8 +314,7 @@ void GeneratorSlowNucleons::GenerateSlow(Int_t charge, Double_t T,
 
 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;
+  Double_t ekin = TMath::Sqrt(p * p + m * m) - m;
   return (p * p * exp(-ekin / T));
 }
 

GeneratorSlowNucleons/GeneratorSlowNucleons.h

@@ -43,11 +43,11 @@ class GeneratorSlowNucleons : public TGenerator {
   //    {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) {
+  virtual void SetTemperature(Float_t t1 = 0.04, Float_t t2 = 0.004) {
     fTemperatureG = t1;
     fTemperatureB = t2;
   }
-  virtual void SetBetaSource(Double_t b1 = 0.05, Double_t b2 = 0.) {
+  virtual void SetBetaSource(Float_t b1 = 0.05, Float_t b2 = 0.) {
     fBetaSourceG = b1;
     fBetaSourceB = b2;
   }

GeneratorSlowNucleons/test.C

@@ -2,19 +2,24 @@ void test()
 {
     GeneratorSlowNucleons* slow = new GeneratorSlowNucleons();
     slow->SetSlowNucleonModel(new SlowNucleonModelExp());
-    slow->SetNcoll(7);
     slow->SetNominalCmsEnergy(13000);
     slow->SetTarget();
+    slow->SetTemperature();
+    slow->SetBetaSource();
+    slow->SetProtonDirection();
+    slow->SetDebug();
+    slow->SetPmax();
     slow->Init();
+
+    slow->SetNcoll(7);
     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());
-            
+            printf(" PDG:%d pt=%1.2f GeV pz=%1.2f GeV\n", p.GetPdgCode(), p.Pt(), p.Pz());
     }
 }

GeneratorSpectators/test.C

@@ -0,0 +1,62 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include <TClonesArray.h>
+#include <TDatabasePDG.h>
+#include <TLorentzVector.h>
+#include <TMath.h>
+#include <TPDGCode.h>
+#include <TParticle.h>
+#include <TParticlePDG.h>
+#include <TRandom.h>
+#include <TH1F.h>
+#include <TCanvas.h>
+#include "GeneratorSpectators.h"
+
+#endif
+
+void test(int nev=1)
+{
+    GeneratorSpectators* spgen = new GeneratorSpectators();
+    spgen->SetParticle(2112); //Neutrons
+    spgen->SetMomentum(2510.);
+    spgen->SetDirection(0, 0., 0., -1.);
+    spgen->SetFermi();
+    spgen->SetDivergence(0.000032); //beam divergence in murad
+    spgen->SetCrossing(0.000060, 2);//beam crossing angle in murad
+
+    spgen->Init();
+
+    double upper=4000.;
+    if(spgen->GetZDirection()<0.) upper = -4000.;
+    TH1F *hpz = new TH1F("hpz", "Spectators p_{z}", 200, 0., upper);
+    TH1F *hpt = new TH1F("hpt", "Spectators p_{t}", 100, 0., 2.);
+
+    for (Int_t iev = 0; iev < nev; iev++) {
+      spgen->GenerateEvent();
+
+      //printf(" \t Importing generated particles \n");
+      auto particles = new TClonesArray("TParticle", 10);
+      auto npart = spgen->ImportParticles(particles, "all");
+      //printf(" - Number of generated particles %5d \n", npart);
+
+      for (Int_t i = 0; i < npart; i++) {
+        auto p = *(TParticle*) (particles->At(i));
+        //printf("(%d) PDG:%5d pt=%1.2f GeV pz=%1.2f GeV\n", iev, p.GetPdgCode(), p.Pt(), p.Pz());
+        hpt->Fill(p.Pt());
+        hpz->Fill(p.Pz());
+      }
+    }
+
+      TCanvas *c = new TCanvas("c","Momentum",0,0,1200,600);
+      c->Divide(2,1);
+      c->cd(1);
+      hpt->SetLineColor(kAzure+10);
+      hpt->SetLineWidth(2);
+      hpt->Draw("e");
+      hpt->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+      c->cd(2);
+      hpz->SetLineColor(kPink);
+      hpz->SetLineWidth(2);
+      hpz->Draw("e");
+      hpz->GetXaxis()->SetTitle("p_{Z} (GeV/c)");
+}