+/**************************************************************************
+ * 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$ */
+
+
+
+//-------------------------------------------------------------
+// Generate Beam-Gas-Events
+// Default: underlying event: pO @ 7 TeV (fixed target)
+//
+// underlying event can be changed by adding generators
+// in the same way as to AliGenCocktail before calling
+// the Init()
+// Author: Jochen Klein
+//-------------------------------------------------------------
+
+#include "AliGenBeamGasNew.h"
+#include "AliGenCocktail.h"
+#include "AliGenCocktailEntry.h"
+#include "AliGenEventHeader.h"
+#include "AliGenCocktailEventHeader.h"
+#include "AliGenHijing.h"
+#include <TList.h>
+#include <TObjArray.h>
+#include "AliCollisionGeometry.h"
+#include "AliRun.h"
+#include "AliMC.h"
+#include "AliStack.h"
+#include "AliHeader.h"
+#include "TParticle.h"
+
+ClassImp(AliGenBeamGasNew)
+
+Float_t EtaToTheta(Float_t arg) { return (180./TMath::Pi())*2.*atan(exp(-arg)); }
+
+AliGenBeamGasNew::AliGenBeamGasNew() :
+ AliGenCocktail(),
+ fItime(0),
+ fTwindow(88e-6),
+ fZwindow(2000),
+ fRate(1e3)
+{
+ printf("AliGenBeamGasNew instatiated!\n");
+}
+
+AliGenBeamGasNew::~AliGenBeamGasNew()
+{
+ delete fEntries;
+ fEntries = 0;
+ fHeader = 0;
+}
+
+void AliGenBeamGasNew::SetTimeWindow(Float_t twindow) { fTwindow = twindow; }
+
+bool AliGenBeamGasNew::SetRate(Float_t rate) {
+ if (rate >= 0) {
+ fRate = rate;
+ return true;
+ } else
+ return false;
+}
+
+void AliGenBeamGasNew::SetZWindow(Float_t zwindow) { fZwindow = zwindow; }
+
+void AliGenBeamGasNew::Init()
+{
+ // Initialisation of the class
+ // if no generators were added before calling Init()
+ // AliGenHijing is added configured for pO
+
+ fVertexSmear = kPerEvent;
+ fVertexSource = kInternal;
+ fRandom = kTRUE;
+
+ // Adding default underlying event in case none was specified
+ // p-O-collision at 7 TeV (fixed target)
+ if (!fEntries) {
+ AliGenHijing *gen = new AliGenHijing(-1);
+ gen->SetEnergyCMS(7000);
+ gen->SetReferenceFrame("LAB");
+ gen->SetProjectile("P",1,1);
+ gen->SetTarget("A",16,8);
+ gen->SetPhiRange(0,360);
+ Float_t thmin = EtaToTheta(8);
+ Float_t thmax = EtaToTheta(-8);
+ gen->SetThetaRange(thmin, thmax);
+ gen->SetOrigin(0,0,0);
+ gen->SetSigma(0,0,0);
+ gen->KeepFullEvent();
+ gen->SetShadowing(1);
+ gen->SetDecaysOff(1);
+ gen->SetSpectators(1);
+ gen->SetSelectAll(1);
+ gen->SetRandomPz(kFALSE);
+ gen->SetPileUpTimeWindow(0.);
+ AddGenerator(gen,"Hijing pO",1);
+ }
+
+ TIter next(fEntries);
+ AliGenCocktailEntry *entry;
+
+ // Loop over generators and initialize
+ while((entry = (AliGenCocktailEntry*)next())) {
+ if (fStack) entry->Generator()->SetStack(fStack);
+ entry->Generator()->Init();
+ }
+
+ next.Reset();
+
+ if (fRandom) {
+ fProb.Set(fNGenerators);
+ next.Reset();
+ Float_t sum = 0.;
+ while((entry = (AliGenCocktailEntry*)next())) {
+ sum += entry->Rate();
+ }
+
+ next.Reset();
+ Int_t i = 0;
+ Float_t psum = 0.;
+ while((entry = (AliGenCocktailEntry*)next())) {
+ psum += entry->Rate() / sum;
+ fProb[i++] = psum;
+ }
+ }
+ next.Reset();
+}
+
+void AliGenBeamGasNew::VertexInternal()
+{
+ // calculation of the interaction vertex for the beam gas event
+ // both spatial and time coordinate are adjusted.
+
+ Float_t random[2];
+ Float_t nbunch;
+ Float_t bunch;
+
+ gRandom->RndmArray(2,random);
+ fVertex[0] = fVertex[1] = 0;
+ fVertex[2] = random[0] * 4000 - 2000;
+ nbunch = fTwindow/25e-9;
+ bunch = floor(2*nbunch * random[1] - nbunch);
+ fItime = fVertex[2] / 100 / 3e8 + bunch * 25e-9;
+}
+
+void AliGenBeamGasNew::Generate()
+{
+//
+// Generate the collisions for one event
+//
+ Int_t numbg = gRandom->Poisson(fRate*fZwindow/100*2*fTwindow);
+ if (numbg == 0) return;
+
+ TIter next(fEntries);
+ AliGenCocktailEntry *entry = 0;
+ AliGenCocktailEntry *preventry = 0;
+ AliGenerator* gen = 0;
+ TLorentzVector v;
+ TArrayF eventVertex;
+ Int_t lastpart=0;
+ Float_t sign;
+
+ if (fHeader) delete fHeader;
+ fHeader = new AliGenCocktailEventHeader("Beamgas Header");
+
+ TObjArray *partArray = gAlice->GetMCApp()->Particles();
+ AliStack *stack = gAlice->GetRunLoader()->Stack();
+
+ for (Int_t l = 0; l < numbg; l++) {
+ Vertex();
+ sign = (gRandom->Rndm() < 0.5)? -1. : 1.;
+ eventVertex.Set(3);
+ for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
+
+ if (!fRandom) {
+ Int_t igen=0;
+ while((entry = (AliGenCocktailEntry*)next())) {
+ if (fUsePerEventRate && (gRandom->Rndm() > entry->Rate())) continue;
+
+ igen++;
+ if (igen ==1) {
+ entry->SetFirst(lastpart);
+ } else {
+ entry->SetFirst((partArray->GetEntriesFast())+1);
+ }
+ //
+ // Handle case in which current generator needs collision geometry from previous generator
+ //
+ gen = entry->Generator();
+ if (gen->NeedsCollisionGeometry())
+ {
+ if (preventry && preventry->Generator()->ProvidesCollisionGeometry())
+ {
+ gen->SetCollisionGeometry(preventry->Generator()->CollisionGeometry());
+ } else {
+ Fatal("Generate()", "No Collision Geometry Provided");
+ }
+ }
+ entry->Generator()->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
+ entry->Generator()->Generate();
+ entry->SetLast(partArray->GetEntriesFast());
+ preventry = entry;
+ }
+ } else {
+ //
+ // Select a generator randomly
+ //
+ Int_t i;
+ Float_t p0 = gRandom->Rndm();
+
+ for (i = 0; i < fNGenerators; i++) {
+ if (p0 < fProb[i]) break;
+ }
+
+ entry = (AliGenCocktailEntry*) fEntries->At(i);
+ entry->SetFirst(lastpart);
+ gen = entry->Generator();
+ entry->Generator()->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
+ entry->Generator()->Generate();
+ entry->SetLast(partArray->GetEntriesFast());
+ }
+
+ for (int k = lastpart; k < stack->GetNprimary(); k++) {
+ stack->Particle(k)->ProductionVertex(v);
+ v[2] *= sign;
+ v[3] = fItime; // ??? +=
+ stack->Particle(k)->SetProductionVertex(v);
+ stack->Particle(k)->Momentum(v);
+ v[2] *= sign;
+ stack->Particle(k)->SetMomentum(v);
+ }
+
+ lastpart = stack->GetNprimary();
+ next.Reset();
+ }
+ // Event Vertex
+ fHeader->SetPrimaryVertex(eventVertex);
+ // fHeader->SetNvertex(numbg);
+ if (fContainer) {
+ fContainer->AddHeader(fHeader);
+ } else {
+ gAlice->SetGenEventHeader(fHeader);
+ }
+}