EVGEN/AliGenBeamGasNew.cxx

   1 //-------------------------------------------------------------
   2 //  Generate Beam-Gas-Events
   3 //  Default: underlying event: pO @ 7 TeV (fixed target)
   4 //
   5 //  underlying event can be changed by adding generators
   6 //  in the same as to AliGenCocktail before calling the Init()
   7 //-------------------------------------------------------------
   8
   9 /* $Id$ */
  10
  11 #include <TList.h>
  12 #include <TObjArray.h>
  13 #include "TParticle.h"
  14
  15 #include "AliGenBeamGasNew.h"
  16 #include "AliGenCocktail.h"
  17 #include "AliGenCocktailEntry.h"
  18 #include "AliGenCocktailEventHeader.h"
  19 #include "../THijing/AliGenHijing.h"
  20 #include "AliCollisionGeometry.h"
  21 #include "AliRun.h"
  22 #include "AliMC.h"
  23 #include "AliStack.h"
  24 #include "AliLog.h"
  25
  26 ClassImp(AliGenBeamGasNew)
  27
  28 Float_t EtaToTheta(Float_t arg) { return (180./TMath::Pi())*2.*atan(exp(-arg)); }
  29
  30 AliGenBeamGasNew::AliGenBeamGasNew() :
  31     AliGenCocktail(),
  32     fTwindow(88.e-6)
  33 {
  34     // Default constructor
  35 }
  36
  37 AliGenBeamGasNew::~AliGenBeamGasNew()
  38 {
  39     // Destructor
  40 }
  41
  42 AliGenBeamGasNew::AliGenBeamGasNew(const AliGenBeamGasNew& /*rhs*/) : AliGenCocktail()
  43 {
  44   AliFatal("Copy Constructor not yet implemented!");
  45 }
  46
  47
  48 void AliGenBeamGasNew::Init()
  49 {
  50   //  Initialisation of the class
  51   //  if no generators were added before calling Init()
  52   //  AliGenHijing is added configured for pO
  53
  54   fVertexSmear = kPerEvent;
  55   fVertexSource = kInternal;
  56   fRandom = kTRUE;
  57
  58   // Adding default underlying event in case none was specified
  59   // p-O-collision at 7 TeV (fixed target)
  60   if (!fEntries) {
  61     AliGenHijing *gen = new AliGenHijing(-1);
  62     gen->SetEnergyCMS(7000);
  63     gen->SetReferenceFrame("LAB");
  64     gen->SetProjectile("P",1,1);
  65     gen->SetTarget("A",16,8);
  66     gen->SetPhiRange(0,360);
  67     Float_t thmin = EtaToTheta(8);
  68     Float_t thmax = EtaToTheta(-8);
  69     gen->SetThetaRange(thmin, thmax);
  70     gen->SetOrigin(0,0,0);
  71     gen->SetSigma(0,0,0);
  72     gen->KeepFullEvent();
  73     gen->SetShadowing(1);
  74     gen->SetDecaysOff(1);
  75     gen->SetSpectators(1);
  76     gen->SetSelectAll(1);
  77     gen->SetRandomPz(kTRUE);
  78     gen->SetPileUpTimeWindow(0.);
  79     AddGenerator(gen,"Hijing pO",1);
  80   }
  81
  82   TIter next(fEntries);
  83   AliGenCocktailEntry *entry;
  84
  85   // Loop over generators and initialize
  86   while((entry = (AliGenCocktailEntry*)next())) {
  87       if (fStack)  entry->Generator()->SetStack(fStack);
  88       entry->Generator()->Init();
  89   }
  90
  91   next.Reset();
  92
  93   if (fRandom) {
  94     fProb.Set(fNGenerators);
  95     next.Reset();
  96     Float_t sum = 0.;
  97     while((entry = (AliGenCocktailEntry*)next())) {
  98       sum += entry->Rate();
  99     }
 100
 101     next.Reset();
 102     Int_t i = 0;
 103     Float_t psum = 0.;
 104     while((entry = (AliGenCocktailEntry*)next())) {
 105       psum +=  entry->Rate() / sum;
 106       fProb[i++] = psum;
 107     }
 108   }
 109   next.Reset();
 110 }
 111
 112 void AliGenBeamGasNew::VertexInternal()
 113 {
 114   // calculation of the interaction vertex for the beam gas event
 115   // both spatial and time coordinate are adjusted.
 116
 117     Float_t random[2];
 118     Float_t nbunch;
 119     Int_t bunch;
 120
 121     Rndm(random,2);
 122     fVertex[0] = fVertex[1] = 0.;
 123     fVertex[2] = random[1] * 4000. - 2000.;
 124     AliInfo(Form("z-Koord.: %13.3f \n", fVertex[2]));
 125     nbunch = 2 * fTwindow / 25.e-9;
 126     bunch = (Int_t) (nbunch * (random[0] - 0.5));
 127     fItime = fVertex[2] / 3.e8 + bunch * 25.e-9;
 128     AliInfo(Form("fItime: %13.3e \n", fItime));
 129 }
 130
 131 void AliGenBeamGasNew::Generate()
 132 {
 133   TIter next(fEntries);
 134   AliGenCocktailEntry *entry = 0;
 135   AliGenCocktailEntry *preventry = 0;
 136   AliGenerator* gen = 0;
 137
 138   if (fHeader) delete fHeader;
 139   fHeader = new AliGenCocktailEventHeader("Beamgas Header");
 140
 141   TObjArray *partArray = gAlice->GetMCApp()->Particles();
 142
 143   TArrayF eventVertex;
 144   int numbg = 0; // variable for counting how many bg interactions happened
 145
 146   do {
 147     numbg++;
 148     Vertex();
 149     eventVertex.Set(3);
 150     for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
 151
 152     if (!fRandom) {
 153       Int_t igen=0;
 154       while((entry = (AliGenCocktailEntry*)next())) {
 155         if (fUsePerEventRate && (gRandom->Rndm() > entry->Rate())) continue;
 156
 157         igen++;
 158         if (igen ==1) {
 159           entry->SetFirst(0);
 160         } else {
 161           entry->SetFirst((partArray->GetEntriesFast())+1);
 162         }
 163         //
 164         //      Handle case in which current generator needs collision geometry from previous generator
 165         //
 166         gen = entry->Generator();
 167         if (gen->NeedsCollisionGeometry())
 168           {
 169             if (preventry && preventry->Generator()->ProvidesCollisionGeometry())
 170               {
 171                 gen->SetCollisionGeometry(preventry->Generator()->CollisionGeometry());
 172               } else {
 173                 Fatal("Generate()", "No Collision Geometry Provided");
 174               }
 175           }
 176         entry->Generator()->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
 177         entry->Generator()->Generate();
 178         entry->SetLast(partArray->GetEntriesFast());
 179         preventry = entry;
 180       }
 181     } else {
 182       //
 183       // Select a generator randomly
 184       //
 185       Int_t i;
 186       Float_t p0 =  gRandom->Rndm();
 187
 188       for (i = 0; i < fNGenerators; i++) {
 189         if (p0 < fProb[i]) break;
 190       }
 191
 192       entry = (AliGenCocktailEntry*) fEntries->At(i);
 193       entry->SetFirst(0);
 194       gen = entry->Generator();
 195       entry->Generator()->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
 196       entry->Generator()->Generate();
 197       entry->SetLast(partArray->GetEntriesFast());
 198     }
 199
 200     AliStack *stack = gAlice->GetRunLoader()->Stack();
 201     TLorentzVector v;
 202     for (int k = 0; k < stack->GetNprimary(); k++) {
 203       stack->Particle(k)->ProductionVertex(v);
 204       v[3] += fItime;
 205       stack->Particle(k)->SetProductionVertex(v);
 206     }
 207
 208     next.Reset();
 209   } while (gRandom->Rndm() < 1. / 5.); // must be adjusted to the rate of bg interactions
 210   // Event Vertex
 211   fHeader->SetPrimaryVertex(eventVertex);
 212   // Total number of produced an stored particles
 213   fHeader->CalcNProduced();
 214
 215   gAlice->SetGenEventHeader(fHeader);
 216 }