THijing/AliGenBeamGasNew.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18
  19
  20 //-------------------------------------------------------------
  21 //  Generate Beam-Gas-Events
  22 //  Default: underlying event: pO @ 7 TeV (fixed target)
  23 //
  24 //  underlying event can be changed by adding generators
  25 //  in the same way as to AliGenCocktail before calling
  26 //  the Init()
  27 //  Author: Jochen Klein
  28 //-------------------------------------------------------------
  29
  30 #include <TList.h>
  31 #include <TObjArray.h>
  32 #include <TParticle.h>
  33
  34 #include "AliGenBeamGasNew.h"
  35 #include "AliGenCocktail.h"
  36 #include "AliGenCocktailEntry.h"
  37 #include "AliGenCocktailEventHeader.h"
  38 #include "../THijing/AliGenHijing.h"
  39 #include "AliCollisionGeometry.h"
  40 #include "AliRun.h"
  41 #include "AliMC.h"
  42 #include "AliStack.h"
  43 #include "AliLog.h"
  44
  45 #include "AliHeader.h"
  46 #include "AliGenEventHeader.h"
  47
  48 ClassImp(AliGenBeamGasNew)
  49
  50 Float_t EtaToTheta(Float_t arg) { return (180./TMath::Pi())*2.*atan(exp(-arg)); }
  51
  52 AliGenBeamGasNew::AliGenBeamGasNew() :
  53     AliGenCocktail(),
  54     fItime(0),
  55     fTwindow(88e-6),
  56     fZwindow(2000),
  57     fRate(1e3)
  58 {
  59   // Default constructor
  60 }
  61
  62 AliGenBeamGasNew::~AliGenBeamGasNew()
  63 {
  64 }
  65
  66 void AliGenBeamGasNew::SetTimeWindow(Float_t twindow) { fTwindow = twindow; }
  67
  68 bool AliGenBeamGasNew::SetRate(Float_t rate) {
  69 // Set the Rate
  70   if (rate >= 0) {
  71     fRate = rate;
  72     return true;
  73   } else
  74     return false;
  75 }
  76
  77 void AliGenBeamGasNew::SetZWindow(Float_t zwindow) { fZwindow = zwindow; }
  78
  79 void AliGenBeamGasNew::Init()
  80 {
  81   //  Initialisation of the class
  82   //  if no generators were added before calling Init()
  83   //  AliGenHijing is added configured for pO
  84
  85   fVertexSmear = kPerEvent;
  86   fVertexSource = kInternal;
  87   fRandom = kTRUE;
  88
  89   // Adding default underlying event in case none was specified
  90   // p-O-collision at 7 TeV (fixed target)
  91   if (!fEntries) {
  92     AliGenHijing *gen = new AliGenHijing(-1);
  93     gen->SetReferenceFrame("LAB");
  94     gen->SetEnergyCMS(7000);
  95     gen->SetProjectile("P",1,1);
  96     gen->SetTarget("A",16,8);
  97     gen->SetPhiRange(0,360);
  98     Float_t thmin = EtaToTheta(8);
  99     Float_t thmax = EtaToTheta(-8);
 100     gen->SetThetaRange(thmin, thmax);
 101     gen->SetOrigin(0,0,0);
 102     gen->SetSigma(0,0,0);
 103     gen->KeepFullEvent();
 104     gen->SetShadowing(1);
 105     gen->SetDecaysOff(1);
 106     gen->SetSpectators(1);
 107     gen->SetSelectAll(1);
 108     gen->SetRandomPz(kFALSE);
 109     gen->SetPileUpTimeWindow(0.);
 110     AddGenerator(gen,"Hijing pO",1);
 111   }
 112
 113   TIter next(fEntries);
 114   AliGenCocktailEntry *entry;
 115
 116   // Loop over generators and initialize
 117   while((entry = (AliGenCocktailEntry*)next())) {
 118     if (fStack)  entry->Generator()->SetStack(fStack);
 119     entry->Generator()->Init();
 120   }
 121
 122   next.Reset();
 123
 124   if (fRandom) {
 125     fProb.Set(fNGenerators);
 126     next.Reset();
 127     Float_t sum = 0.;
 128     while((entry = (AliGenCocktailEntry*)next())) {
 129       sum += entry->Rate();
 130     }
 131
 132     next.Reset();
 133     Int_t i = 0;
 134     Float_t psum = 0.;
 135     while((entry = (AliGenCocktailEntry*)next())) {
 136       psum +=  entry->Rate() / sum;
 137       fProb[i++] = psum;
 138     }
 139   }
 140   next.Reset();
 141 }
 142
 143 void AliGenBeamGasNew::VertexInternal()
 144 {
 145   // calculation of the interaction vertex for the beam gas event
 146   // both spatial and time coordinate are adjusted.
 147
 148   Float_t random[2];
 149   Float_t nbunch;
 150   Float_t bunch;
 151
 152   gRandom->RndmArray(2,random);
 153   fVertex[0] = fVertex[1] = 0;
 154   fVertex[2] = random[0] * 2 * fZwindow - fZwindow;
 155   nbunch = fTwindow/25e-9;
 156   bunch = floor(2*nbunch * random[1] - nbunch);
 157   fItime = fVertex[2] / 100 / 3e8 + bunch * 25e-9;
 158   AliInfo(Form("fItime: %13.3e \n", fItime));
 159 }
 160
 161 void AliGenBeamGasNew::Generate()
 162 {
 163 //
 164 // Generate the collisions for one event
 165 //
 166   Int_t numbg = gRandom->Poisson(fRate*fZwindow/100*2*fTwindow);
 167   if (numbg == 0) return;
 168
 169   TIter next(fEntries);
 170   AliGenCocktailEntry *entry = 0;
 171   AliGenCocktailEntry *preventry = 0;
 172   AliGenerator* gen = 0;
 173   TLorentzVector v;
 174   TArrayF eventVertex;
 175   Int_t lastpart=0;
 176   Float_t sign;
 177
 178   if (fHeader) delete fHeader;
 179   fHeader = new AliGenCocktailEventHeader("Beamgas Header");
 180
 181   const TObjArray *partArray = gAlice->GetMCApp()->Particles();
 182   AliStack *stack = AliRunLoader::Instance()->Stack();
 183
 184   for (Int_t l = 0; l < numbg; l++) {
 185     Vertex();
 186     sign = (gRandom->Rndm() < 0.5)? -1. : 1.;
 187     eventVertex.Set(3);
 188     for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
 189
 190     if (!fRandom) {
 191       Int_t igen=0;
 192       while((entry = (AliGenCocktailEntry*)next())) {
 193         if (fUsePerEventRate && (gRandom->Rndm() > entry->Rate())) continue;
 194
 195         igen++;
 196         if (igen ==1) {
 197           entry->SetFirst(lastpart);
 198         } else {
 199           entry->SetFirst((partArray->GetEntriesFast())+1);
 200         }
 201         //
 202         //      Handle case in which current generator needs collision geometry from previous generator
 203         //
 204         gen = entry->Generator();
 205         if (gen->NeedsCollisionGeometry())
 206           {
 207             if (preventry && preventry->Generator()->ProvidesCollisionGeometry())
 208               {
 209                 gen->SetCollisionGeometry(preventry->Generator()->CollisionGeometry());
 210               } else {
 211                 Fatal("Generate()", "No Collision Geometry Provided");
 212               }
 213           }
 214         entry->Generator()->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
 215         entry->Generator()->Generate();
 216         entry->SetLast(partArray->GetEntriesFast());
 217         preventry = entry;
 218       }
 219     } else {
 220       //
 221       // Select a generator randomly
 222       //
 223       Int_t i = 0;
 224       Float_t p0 =  gRandom->Rndm();
 225
 226       for (i = 0; i < fNGenerators; i++) {
 227         if (p0 < fProb[i]) break;
 228       }
 229
 230       entry = (AliGenCocktailEntry*) fEntries->At(i);
 231       entry->SetFirst(lastpart);
 232       gen = entry->Generator();
 233       gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
 234       gen->Generate();
 235       entry->SetLast(partArray->GetEntriesFast());
 236     }
 237
 238     for (Int_t k = lastpart; k < stack->GetNprimary(); k++) {
 239       stack->Particle(k)->ProductionVertex(v);
 240       v[2] *= sign;
 241       eventVertex[2] *= sign;
 242       v[3] = fItime;  // ??? +=
 243       stack->Particle(k)->SetProductionVertex(v);
 244       stack->Particle(k)->Momentum(v);
 245       v[2] *= sign;
 246       stack->Particle(k)->SetMomentum(v);
 247     }
 248
 249     ((AliGenEventHeader*) fHeader->GetHeaders()->Last())->SetPrimaryVertex(eventVertex);
 250
 251     lastpart = stack->GetNprimary();
 252     next.Reset();
 253   }
 254
 255   fHeader->CalcNProduced();
 256
 257   if (fContainer) {
 258     fContainer->AddHeader(fHeader);
 259   } else {
 260     gAlice->SetGenEventHeader(fHeader);
 261   }
 262 }