]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EVGEN/AliGenSlowNucleons.cxx
Script to create a random bad channel map.
[u/mrichter/AliRoot.git] / EVGEN / AliGenSlowNucleons.cxx
index ba0388521cbc45d33b47971a1ef4cbdc7adb8fa2..16cb40b0a71e1c54316af583b88e94e1baa00e66 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*
-$Log$
-*/
+/* $Id$ */
 
-/*
-  Generator for slow nucluons in pA interactions. 
-  Source is modelled by a relativistic Maxwell distributions.
-  Original code by  Ferenc Sikler  <sikler@rmki.kfki.hu>
- */
+//
+//  Generator for slow nucleons in pA interactions. 
+//  Source is modelled by a relativistic Maxwell distributions.
+//  This class cooparates with AliCollisionGeometry if used inside AliGenCocktail.
+//  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  <sikler@rmki.kfki.hu>
+// 
 
 #include <TDatabasePDG.h>
 #include <TPDGCode.h>
 #include <TH2F.h>
+#include <TH1F.h>
+#include <TF1.h>
+#include <TCanvas.h>
 
 #include "AliCollisionGeometry.h"
 #include "AliGenSlowNucleons.h"
 #include "AliSlowNucleonModel.h"
 
- ClassImp(AliGenSlowNucleons)
+ClassImp(AliGenSlowNucleons)
+
     
- AliGenSlowNucleons::AliGenSlowNucleons():AliGenerator(-1)
+AliGenSlowNucleons::AliGenSlowNucleons()
+    :AliGenerator(-1),
+     fCMS(0.),
+     fMomentum(0.),
+     fBeta(0.),
+     fPmax (0.),
+     fATarget (0.),
+     fZTarget (0.),
+     fCharge(0),
+     fProtonDirection(0.),
+     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(),
+     fSlowNucleonModel(0)
 {
 // Default constructor
-    fSlowNucleonModel = 0;
     fCollisionGeometry = 0;
 }
 
 AliGenSlowNucleons::AliGenSlowNucleons(Int_t npart)
-    :AliGenerator(npart)
+    :AliGenerator(npart),
+     fCMS(14000.),
+     fMomentum(0.),
+     fBeta(0.),
+     fPmax (10.),
+     fATarget (208.),
+     fZTarget (82.),
+     fCharge(1),
+     fProtonDirection(0.),
+     fTemperatureG(0.04), 
+     fBetaSourceG(0.05),
+     fTemperatureB(0.004),
+     fBetaSourceB(0.),
+     fNgp(0),
+     fNgn(0),
+     fNbp(0),
+     fNbn(0),
+     fDebug(0),
+     fDebugHist1(0),
+     fDebugHist2(0),
+     fThetaDistribution(),
+     fCosThetaGrayHist(),
+     fCosTheta(),
+     fSlowNucleonModel(new AliSlowNucleonModel())
 {
 // Constructor
     fName  = "SlowNucleons";
     fTitle = "Generator for gray particles in pA collisions";
-    SetPmax();
-    SetTarget();
-    SetNominalCmsEnergy();
-    SetCharge();
-    SetTemperature();
-    SetBetaSource();
-    fSlowNucleonModel = new AliSlowNucleonModel();
     fCollisionGeometry = 0;
-    fDebug = 0;
 }
 
 //____________________________________________________________
@@ -64,6 +107,10 @@ AliGenSlowNucleons::~AliGenSlowNucleons()
     delete  fSlowNucleonModel;
 }
 
+void AliGenSlowNucleons::SetProtonDirection(Float_t dir) {
+// Set direction of the proton to change between pA (1) and Ap (-1)
+  fProtonDirection = dir / TMath::Abs(dir);
+}
 
 void AliGenSlowNucleons::Init()
 {
@@ -74,14 +121,33 @@ void AliGenSlowNucleons::Init()
     fMomentum = fCMS/2. * fZTarget / fATarget;
     fBeta     = fMomentum / TMath::Sqrt(kMass * kMass + fMomentum * fMomentum);
     if (fDebug) {
-       fDebugHist = new TH2F("DebugHist", "N_heavy vs nu", 50, 0., 50., 50, 0., 50.);
+       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.);
     }
 }
 
 void AliGenSlowNucleons::FinishRun()
 {
+// End of run action
+// Show histogram for debugging if requested.
     if (fDebug) {
-       fDebugHist->Draw();
+       TCanvas *c = new TCanvas("c","Canvas 1",400,10,600,700);
+       c->Divide(2,1);
+       c->cd(1);
+       fDebugHist1->Draw();
+       c->cd(2);
+       fDebugHist2->Draw();
+       c->cd(3);
+       fCosThetaGrayHist->Draw();
     }
 }
 
@@ -94,40 +160,43 @@ void AliGenSlowNucleons::Generate()
   //
   // Communication with Gray Particle Model 
   // 
-    Int_t ngp, ngn, nbp, nbn;
-
-    Float_t b = fCollisionGeometry->ImpactParameter();
-    Int_t  nn   = fCollisionGeometry->NN();
-    Int_t  nwn  = fCollisionGeometry->NwN();
-    Int_t  nnw  = fCollisionGeometry->NNw();
-    Int_t  nwnw = fCollisionGeometry->NwNw();
-
-    printf("AliGenSlowNucleons: Impact parameter from Collision Geometry %f %d %d %d %d\n", 
-          b, nn, nwn, nnw, nwnw);
-        
-    fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, ngp, ngn, nbp, nbn);
-
-    if (fDebug) {
-       printf("nucleons %d %d %d %d \n", ngp, ngn, nbp, nbn);
-       
-       fDebugHist->Fill(Float_t(ngp + ngn + nbp + nbn), fCollisionGeometry->NwN(), 1.);
-    }
-    
+    if (fCollisionGeometry) {
+       Float_t b   = fCollisionGeometry->ImpactParameter();
+       Int_t  nn   = fCollisionGeometry->NN();
+       Int_t  nwn  = fCollisionGeometry->NwN();
+       Int_t  nnw  = fCollisionGeometry->NNw();
+       Int_t  nwnw = fCollisionGeometry->NwNw();
+
+       fSlowNucleonModel->GetNumberOfSlowNucleons(fCollisionGeometry, fNgp, fNgn, fNbp, fNbn);
+       if (fDebug) {
+           printf("Nucleons %d %d %d %d \n", fNgp, fNgn, fNbp, fNbn);
+           fDebugHist1->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), fCollisionGeometry->NwN(), 1.);
+           fDebugHist2->Fill(Float_t(fNgp + fNgn + fNbp + fNbn), b, 1.);
+           printf("AliGenSlowNucleons: Impact parameter from Collision Geometry %f %d %d %d %d\n", 
+                  b, nn, nwn, nnw, nwnw);
+       }
+    }     
 
    //
-    Float_t p[3];
+    Float_t p[3], theta=0;
     Float_t origin[3] = {0., 0., 0.};
     Float_t polar [3] = {0., 0., 0.};    
-    Int_t nt, i;
+    Int_t nt, i, j;
     Int_t kf;
+    
+    if(fVertexSmear == kPerEvent) {
+       Vertex();
+       for (j=0; j < 3; j++) origin[j] = fVertex[j];
+    } // if kPerEvent
 //
 //  Gray protons
 //
     fCharge = 1;
     kf = kProton;    
-    for(i = 0; i < ngp; i++) {
-       GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p);
-       SetTrack(fTrackIt, -1, kf, p, origin, polar,
+    for(i = 0; i < fNgp; i++) {
+       GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
+       if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
+       PushTrack(fTrackIt, -1, kf, p, origin, polar,
                 0., kPNoProcess, nt, 1.);
        KeepTrack(nt);
     }
@@ -136,9 +205,10 @@ void AliGenSlowNucleons::Generate()
 //
     fCharge = 0;
     kf = kNeutron;    
-    for(i = 0; i < ngn; i++) {
-       GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p);
-       SetTrack(fTrackIt, -1, kf, p, origin, polar,
+    for(i = 0; i < fNgn; i++) {
+       GenerateSlow(fCharge, fTemperatureG, fBetaSourceG, p, theta);
+       if (fDebug) fCosThetaGrayHist->Fill(TMath::Cos(theta));
+       PushTrack(fTrackIt, -1, kf, p, origin, polar,
                 0., kPNoProcess, nt, 1.);
        KeepTrack(nt);
     }
@@ -147,9 +217,9 @@ void AliGenSlowNucleons::Generate()
 //
     fCharge = 1;
     kf = kProton;    
-    for(i = 0; i < nbp; i++) {
-       GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p);
-       SetTrack(fTrackIt, -1, kf, p, origin, polar,
+    for(i = 0; i < fNbp; i++) {
+       GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
+       PushTrack(fTrackIt, -1, kf, p, origin, polar,
                 0., kPNoProcess, nt, 1.);
        KeepTrack(nt);
     }
@@ -158,9 +228,9 @@ void AliGenSlowNucleons::Generate()
 //
     fCharge = 0;
     kf = kNeutron;    
-    for(i = 0; i < nbn; i++) {
-       GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p);
-       SetTrack(fTrackIt, -1, kf, p, origin, polar,
+    for(i = 0; i < fNbn; i++) {
+       GenerateSlow(fCharge, fTemperatureB, fBetaSourceB, p, theta);
+       PushTrack(fTrackIt, -1, kf, p, origin, polar,
                 0., kPNoProcess, nt, 1.);
        KeepTrack(nt);
     }
@@ -169,16 +239,17 @@ void AliGenSlowNucleons::Generate()
 
 
 
-void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T, Double_t beta, Float_t* q)
+void AliGenSlowNucleons::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]    
 */
 
-{
- Double_t m, pmax, p, f, theta, phi;
+ 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();
@@ -202,8 +273,14 @@ void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T, Double_t beta, F
  }
  while(f < Rndm());
 
- /* Spherical symmetric emission */
- theta = TMath::ACos(2. * Rndm() - 1.);
+ /* Spherical symmetric emission for black particles (beta=0)*/
+ if(beta==0 || fThetaDistribution==0) theta = TMath::ACos(2. * 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() * Rndm();
 
  /* Determine momentum components in system of the moving source */
@@ -217,6 +294,7 @@ void AliGenSlowNucleons::GenerateSlow(Int_t charge, Double_t T, Double_t beta, F
 
  /* Transform to laboratory system */
  Lorentz(m, fBeta, q);
+ q[2] *= fProtonDirection; 
 }
 
 Double_t AliGenSlowNucleons::Maxwell(Double_t m, Double_t p, Double_t T)
@@ -236,10 +314,3 @@ void AliGenSlowNucleons::Lorentz(Double_t m, Double_t beta, Float_t* q)
     Double_t energy = sqrt(m*m + q[0]*q[0] + q[1]*q[1] + q[2]*q[2]);
     q[2] = gamma * (q[2] + beta*energy);
 }
-
-
-
-
-
-
-