New class added.
[u/mrichter/AliRoot.git] / EVGEN / AliGenPileup.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 //-------------------------------------------------------------------------
17 //                          Class AliGenPileup
18 //   This is a generator of beam-beam pileup.
19 //   It generates interactions within 3 orbits (+-1) around
20 //   the trigger event. The trigger event itself is chosen
21 //   randomly among the bunch crossings within the central orbit.
22 //   The user can decide whenever to include in the simulation the
23 //   "trigger" interaction or not. This is handled by the
24 //   GenerateTrigInteraction(Bool_t flag) method.
25 //   In the case the trigger interaction is included, it is
26 //   generated using the same settings (vertex smear for example) as
27 //   the pileup events.
28 //   In case the trigger simulation is not included, the user can make
29 //   a cocktail of generator used to produce the trigger interaction and
30 //   AliGenPileup. In this case in order to avoid a fake increase of the rate around the
31 //   trigger, the number of background events within the bunch
32 //   crossing of the trigger is readuced by one.
33 //   The beam profile (the list of the active bunch crossings) can be
34 //   controlled via the SetBCMask(const char *mask) method. The syntax
35 //   follows the one in AliTriggerBCMask class. For example:
36 //   "3564H" would mean that all the bunch corssings within the orbit
37 //   are aloowed (which is of course unphysical). In case one wants to simulate
38 //   one-bunch-crossing-per-orbit scenario, the way to do it is to put something like:
39 //   "1H3563L" or similar.
40 //   The SetGenerator(AliGenerator *generator, Float_t rate) method is
41 //   used in order to define the generator to be used. The second argument is the pileup
42 //   rate in terms of #_of_interactions/bunch-crossing = sigma_tot * luminosity.
43 //   The pileup generation time window can be set via
44 //   AliGenerator::SetPileUpTimeWindow(Float_t pileUpTimeW) method. By the default the
45 //   window is set to 88micros (= TPC readout window).
46 //      
47 // cvetan.cheshkov@cern.ch  9/12/2008
48 //-------------------------------------------------------------------------
49
50 #include <TParticle.h>
51 #include <TFormula.h>
52
53 #include "AliGenPileup.h"
54 #include "AliLog.h"
55 #include "AliGenCocktailEventHeader.h"
56 #include "AliGenCocktailEntry.h"
57 #include "AliRun.h"
58 #include "AliStack.h"
59
60 ClassImp(AliGenPileup)
61
62 AliGenPileup::AliGenPileup():
63   AliGenCocktail(),
64   fBCMask("bcm","3564H"),
65   fGenTrig(kFALSE),
66   fFlag(kFALSE)
67 {
68 // Constructor
69 // The pileup time window is by default
70 // set to the TPC readout one
71     fName = "Pileup";
72     fTitle= "Beam-beam pileup";
73
74     fPileUpTimeWindow = 88e-6;
75 }
76
77 AliGenPileup::~AliGenPileup()
78 {
79 // Destructor
80 }
81
82 void AliGenPileup::SetGenerator(AliGenerator *generator, Float_t rate, Bool_t flag)
83 {
84   // The method sets the geenrator to be used
85   // for pileup simulation.
86   // The second argument is the pileup rate in terms of
87   // #_of_interactions/bunch-crossing = sigma_tot * luminosity.
88   // There is a protection in case the generator was already set.
89   if (fEntries) {
90     if (FirstGenerator()) {
91       AliError("Pileup generator has been already set! Nothing done");
92       return;
93     }
94   }
95   AddGenerator(generator,"pileup generator",rate);
96   fFlag = flag;
97 }
98
99 void AliGenPileup::AddGenerator(AliGenerator *Generator,
100                                 const char* Name,
101                                 Float_t RateExp , TFormula* /*form*/)
102 {
103   // The method used to add the pileup generator
104   // in the cocktail list.
105   // The method is protected in order to avoid
106   // its misusage
107   AliGenCocktail::AddGenerator(Generator,Name,RateExp);
108 }
109
110 Bool_t AliGenPileup::SetBCMask(const char *mask)
111 {
112   // Set the active bunch-crossings that
113   // will be included in the pileup
114   // simulation. For more details on the
115   // syntax of the mask - see
116   // STEER/AliTriggerBCMask.* and the comments
117   // in the header of this file
118   return fBCMask.SetMask(mask);
119 }
120
121 void AliGenPileup::Generate()
122 {
123   //
124   // Generate pileup event 
125   // For details see the coments inline
126
127   // Check that the pileup generator is correctly set
128   AliGenCocktailEntry *entry = FirstGenerator();
129   if (!entry) {
130     AliFatal("No pileup generator entry is found!");
131   }
132
133   AliGenerator *gen = entry->Generator();
134   if (!gen) {
135     AliFatal("No pileup generator specified!");
136   }
137   else if (gen->NeedsCollisionGeometry()) {
138     AliFatal("No Collision Geometry Provided");
139   }
140
141   // Check that the pileup rate is correctly set
142   Float_t rate = entry->Rate();
143   if (rate <= 0) {
144     AliFatal(Form("Invalid rate value: %f",rate));
145   }
146
147   // Create cocktail header
148   if (fHeader) delete fHeader;
149   fHeader = new AliGenCocktailEventHeader("Pileup Cocktail Header");
150
151   // Generate time of all
152   // the collisions within one orbit
153   Int_t *nIntBC = new Int_t[3*AliTriggerBCMask::kNBits];
154   Int_t *indexBC = new Int_t[3*AliTriggerBCMask::kNBits];
155   Int_t nTotBC = 0;
156   while (nTotBC == 0) {
157     for(Int_t iBC = 0; iBC <  AliTriggerBCMask::kNBits; iBC++) {
158
159       if (!fBCMask.GetMask(iBC)) continue;
160
161       //      Int_t nInteractions = gRandom->Poisson(rate);
162       Int_t nInteractions;
163       if (!fFlag) 
164         nInteractions = gRandom->Poisson(rate);
165       else 
166         nInteractions = TMath::Nint(rate) + 1;
167
168       if (nInteractions == 0) continue;
169
170       nIntBC[nTotBC] = nInteractions;
171       indexBC[nTotBC] = iBC;
172       nTotBC++;
173     }
174   }
175
176   // Select the bunch crossing for triggered event
177   Int_t iTrgBC = gRandom->Integer(nTotBC);
178   // Subtract one from the number of events
179   // generated within this bc (only in case
180   // the user disabled the generation of the trigger
181   // interaction)
182   if (!fGenTrig) nIntBC[iTrgBC]--;
183
184   // Remove bunch crossings outside pileup
185   // time window
186   for(Int_t iBC = 0; iBC <  nTotBC; iBC++) {
187    if (TMath::Abs(25e-9*(indexBC[iBC]-indexBC[iTrgBC])) > fPileUpTimeWindow)
188      nIntBC[iBC] = 0;
189   }
190
191   // Generate the two orbits around the central one
192   // taking into account the pileup time window
193   for(Int_t iBC = 0; iBC <  AliTriggerBCMask::kNBits; iBC++) {
194
195     if (!fBCMask.GetMask(iBC)) continue;
196
197     if (TMath::Abs(25e-9*(iBC-AliTriggerBCMask::kNBits-indexBC[iTrgBC])) > fPileUpTimeWindow) continue;
198
199     Int_t nInteractions = gRandom->Poisson(rate);
200     if (nInteractions == 0) continue;
201
202     nIntBC[nTotBC] = nInteractions;
203     indexBC[nTotBC] = iBC-AliTriggerBCMask::kNBits;
204     nTotBC++;
205   }
206   for(Int_t iBC = 0; iBC <  AliTriggerBCMask::kNBits; iBC++) {
207
208     if (!fBCMask.GetMask(iBC)) continue;
209
210     if (TMath::Abs(25e-9*(iBC+AliTriggerBCMask::kNBits-indexBC[iTrgBC])) > fPileUpTimeWindow) continue;
211
212     Int_t nInteractions = gRandom->Poisson(rate);
213     if (nInteractions == 0) continue;
214
215     nIntBC[nTotBC] = nInteractions;
216     indexBC[nTotBC] = iBC+AliTriggerBCMask::kNBits;
217     nTotBC++;
218   }
219
220   // Loop over the generated collision times, call the generator
221   // and correct the partcile times in the stack
222   AliStack *stack = AliRunLoader::Instance()->Stack();
223   Int_t lastpart=0;
224   entry->SetFirst(lastpart);
225
226   for(Int_t iBC = 0; iBC <  nTotBC; iBC++) {
227     Float_t deltat = 25e-9*(indexBC[iBC] - indexBC[iTrgBC]);
228     for (Int_t i = 0; i < nIntBC[iBC]; i++) {
229       //  Generate the vertex position and time
230       Vertex();
231       TArrayF eventVertex(3);
232       for (Int_t j=0; j < 3; j++) eventVertex[j] = fVertex[j];
233       Double_t vTime = deltat + gRandom->Gaus(0,fOsigma[2]/TMath::Ccgs());
234     
235       gen->SetVertex(fVertex.At(0), fVertex.At(1), fVertex.At(2));
236       gen->Generate();
237
238       for (Int_t k = lastpart; k < stack->GetNprimary(); k++) {
239         TLorentzVector v;
240         stack->Particle(k)->ProductionVertex(v);
241         v[3] = vTime;
242         stack->Particle(k)->SetProductionVertex(v);
243       }
244       lastpart = stack->GetNprimary();
245
246       // Store the interaction header in the container of the headers
247       ((AliGenEventHeader*) fHeader->GetHeaders()->Last())->SetPrimaryVertex(eventVertex);
248       ((AliGenEventHeader*) fHeader->GetHeaders()->Last())->SetInteractionTime(vTime);
249     }
250   }
251   delete [] nIntBC;
252   delete [] indexBC;
253
254   entry->SetLast(stack->GetNprimary());
255
256   fHeader->CalcNProduced();
257
258   if (fContainer) {
259     fContainer->AddHeader(fHeader);
260   } else {
261     gAlice->SetGenEventHeader(fHeader); 
262   }
263  
264 }
265
266 void AliGenPileup::SetRandomise(Bool_t /*flag*/)
267 {
268   // This setting is not implemented in
269   // case of pileup generation
270   // So the method gives an warning and exits
271   AliWarning("This setting has no effect on the generator!");
272 }
273
274 void AliGenPileup::UsePerEventRates()
275 {
276   // This setting is not implemented in
277   // case of pileup generation
278   // So the method gives an warning and exits
279   AliWarning("This setting has no effect on the generator!");
280 }