Merge branch 'MyDevBranch'
[u/mrichter/AliRoot.git] / ACORDE / AliGenACORDE.cxx
CommitLineData
b86e74f5 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// Contain parametrizations to generate atmospheric muons, and also
21// to generate single muons and muon bundles at surface level.
22//
23//Begin_Html
24/*
25<img src="picts/AliGenACORDEClass.gif">
26</pre>
27<br clear=left>
28<font size=+2 color=red>
29<p>The responsible person for this module is
30<a href="mailto:Enrique.Gamez.Flores@cern.ch">Enrique Gamez</a>.
31</font>
32<pre>
33*/
34//End_Html
35//
36/////////////////////////////////////////////////////////////////////////////
37
38#include "AliGenACORDE.h"
39
40#include <TMCProcess.h>
41#include <TPDGCode.h>
42#include <TClonesArray.h>
43#include <TF1.h>
44#include <TH1F.h>
45
46#include "AliRun.h"
47#include "AliConst.h"
48
49ClassImp(AliGenACORDE)
50
51//_____________________________________________________________________________
52AliGenACORDE::AliGenACORDE()
53 : AliGenerator(),
54 fIpart(0),
55 fCRMode(kSingleMuons),
56 fCRModeName(0),
57 fXwidth(0),
58 fNx(1),
59 fZwidth(0),
60 fNz(1),
61 fMuonGrid(kFALSE),
62 fZenithMin(0),
63 fZenithMax(0),
64 fAzimuthMin(0),
65 fAzimuthMax(0),
66 fPRange(0),
67 fPResolution(1),
68 fAp(0),
69 fMomentumDist(0),
70 fUnfoldedMomentumDist(0),
71 fZenithDist(0),
0e95e9da 72 fPDist(0),
73 fNParticles(0)
b86e74f5 74{
75 //
76 // Default ctor.
77 //
78}
79
80//_____________________________________________________________________________
81AliGenACORDE::AliGenACORDE(Int_t npart)
82 : AliGenerator(npart),
83 fIpart(kMuonMinus),
84 fCRMode(kSingleMuons),
85 fCRModeName(0),
86 fXwidth(0),
87 fNx(1),
88 fZwidth(0),
89 fNz(1),
90 fMuonGrid(kFALSE),
91 fZenithMin(0),
92 fZenithMax(0),
93 fAzimuthMin(0),
94 fAzimuthMax(0),
95 fPRange(0),
96 fPResolution(1),
97 fAp(0),
98 fMomentumDist(0),
99 fUnfoldedMomentumDist(0),
100 fZenithDist(0),
1881b46e 101 fPDist(0),
102 fNParticles(0)
b86e74f5 103{
104 //
105 // Standard ctor.
106 //
107 fName = "ACORDE";
108 fTitle = "Cosmic Muons generator";
109
110 // Set the origin above the vertex, on the surface.
111 fOrigin[0] = 0.;
112 fOrigin[1] = AliACORDEConstants::Instance()->Depth(); // At the surface by default.
113 fOrigin[2] = 0.;
114}
115
116//_____________________________________________________________________________
b86e74f5 117AliGenACORDE::~AliGenACORDE()
118{
119 //
120 // Default dtor.
121 //
122 if ( fPDist ) {fPDist->Delete(); delete fPDist; fPDist = 0;}
123 if ( fUnfoldedMomentumDist ) delete fUnfoldedMomentumDist;
124 if ( fMomentumDist ) delete fMomentumDist;
125 if ( fAp ) delete fAp;
126 if ( fCRModeName ) delete fCRModeName;
127}
128
129//_____________________________________________________________________________
130void AliGenACORDE::Generate()
131{
132 //
133 // Generate on one trigger
134 // Call the respective method inside the loop for the number
135 // of tracks per trigger.
136
1881b46e 137 for (Int_t i = 0; i < fNParticles; i++ ) {
b86e74f5 138
139 if ( fCRMode == kMuonBundle ) {
140 this->GenerateOneMuonBundle();
141
142 } else if ( fCRMode == kSingleMuons ) {
143 this->GenerateOneSingleMuon(kTRUE);
144
145 } else {
146 // Generate only single muons following the parametrizations
147 // for atmospheric muons.
148 this->GenerateOneSingleMuon();
149
150 }
151
152 }
153}
154
155//_____________________________________________________________________________
156void AliGenACORDE::Init()
157{
158 //
159 // Initialize some internal methods.
160 //
161
1881b46e 162
163 printf("**************************************************************\n");
164 printf("<<< *** Starting the AliGenACORDE cosmic generator ******** >>>\n");
165 printf("<<< *** No. of muons generated at the surface of P2: %d, * >>>\n",fNParticles);
166 printf("**************************************************************\n");
167
b86e74f5 168 // Determine some specific data members.
169 fPRange = TMath::Abs(fPMax-fPMin);
170
171 if ( fCRMode == kSingleMuons ) {
172 fCRModeName = new TString("Single Muons");
173 // Initialisation, check consistency of selected ranges
174 if(TestBit(kPtRange)&&TestBit(kMomentumRange))
175 Fatal("Init","You should not set the momentum range and the pt range!");
176
177 if((!TestBit(kPtRange))&&(!TestBit(kMomentumRange)))
178 Fatal("Init","You should set either the momentum or the pt range!");
179
180 } else if ( fCRMode == kMuonBundle ) {
181 fCRModeName = new TString("Muon Bundles");
182
183 } else if ( fCRMode == kMuonFlux ) {
184 fCRModeName = new TString("Muon Fluxes");
185 // Initialize the ditribution functions.
186 this->InitMomentumGeneration();
187 this->InitZenithalAngleGeneration();
188
189 } else {
190 Fatal("Generate", "Generation Mode unknown!\n");
191
192 }
193
194}
195
196//____________________________________________________________________________
197void AliGenACORDE::GenerateOneSingleMuon(Bool_t withFlatMomentum)
198{
199 //
200 // Generate One Single Muon
201 // This method will generate only one muon.
202 // The momentum will be randomly flat distributed if
203 // the paremeter "withFlatMomentum" is set to kTRUE,
204 // otherwise the momentum will generate acordingly the parametrization
205 // given by
206 // and adpted from Bruno Alessandro's implementation with the
207 // CERNLIB to AliRoot.
208 // The "withFlatMomentum" parameter also will be used to generate
209 // the muons with a flat Zenithal angle distribution.
210 // Do the smearing here, so that means per track.
211
212 Float_t polar[3]= {0,0,0}; // Polarization parameters
213 Float_t origin[3];
214 Int_t nt;
215 Float_t p[3];
216 Float_t pmom, pt;
217 Float_t random[6];
218
219 // Take the azimuth random.
220 Rndm(random, 2);
221 Float_t azimuth = fAzimuthMin + (fAzimuthMax-fAzimuthMin)*random[0];
222 Float_t zenith = fZenithMin + (fZenithMax - fZenithMin)*random[1];
223
224 if ( withFlatMomentum ) {
225 Rndm(random,3);
226 if(TestBit(kMomentumRange)) {
227 pmom = -( fPMin + random[0]*(fPMax - fPMin) ); // always downwards.
228 pt = pmom*TMath::Sin(zenith*kDegrad);
229 } else {
230 pt = -( fPtMin + random[1]*(fPtMax - fPtMin)); // always downwards.
231 pmom = pt/TMath::Sin(zenith*kDegrad);
232 }
233
234 } else {
235 if ( fMomentumDist ) {
236 pmom = -this->GetMomentum(); // Always downwards.
237 } else {
238 pmom = -fPMin;
239 }
240 zenith = this->GetZenithAngle(pmom); // In degrees
241 pt = pmom*TMath::Sin(zenith*kDegrad);
242 }
243
244 p[0] = pt*TMath::Sin(azimuth*kDegrad);
245 p[1] = pmom*TMath::Cos(zenith*kDegrad);
246 p[2] = pt*TMath::Cos(azimuth*kDegrad);
247
248 // Finaly the origin, with the smearing
249 Rndm(random,6);
250 origin[0] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
251 TMath::Sin(azimuth*kDegrad)
252 + fOsigma[0]* TMath::Cos(2*random[0]*TMath::Pi())*TMath::Sqrt(-2*TMath::Log(random[1]));
253
254 origin[1] = AliACORDEConstants::Instance()->Depth();
255
256 origin[2] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
257 TMath::Cos(azimuth*kDegrad)
258 + fOsigma[2]* TMath::Cos(2*random[2]*TMath::Pi())*TMath::Sqrt(-2*TMath::Log(random[3]));
259
260 // Put the track on the stack.
261 PushTrack(fTrackIt,-1,fIpart,p,origin,polar,0,kPPrimary,nt);
262
263}
264
265//____________________________________________________________________________
266void AliGenACORDE::GenerateOneMuonBundle()
267{
268 //
269 // Generate One Muon Bundle method
270 // This method will generate a bunch of muons following the
271 // procedure of the AliGenScan class.
272 // These muons will be generated with flat momentum.
273
274 Float_t polar[3]= {0,0,0}; // Polarization parameters
275 Float_t origin[3];
276 Float_t p[3];
277 Int_t nt;
278 Float_t pmom;
279 Float_t random[6];
280
281 Rndm(random, 3);
282 Float_t zenith = fZenithMin + (fZenithMax - fZenithMin)*random[1];
283 Float_t azimuth = fAzimuthMin + (fAzimuthMax-fAzimuthMin)*random[2];
284 //Float_t zenith = 10;
285 //Float_t azimuth = 30;
286
287 // Generate the kinematics a la AliGenScan (Andreas Morchs)
288 Float_t dx, dz;
289 if ( fNx > 0 ) {
290 dx = fXwidth/fNx;
291 } else {
292 dx = 1e10;
293 //dx = 100.;
294 }
295
296 if ( fNz > 0 ) {
297 dz = fZwidth/fNz;
298 } else {
299 dz = 1e10;
300 //dz = 100.;
301 }
302
303 origin[0] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
304 TMath::Sin(azimuth*kDegrad);
305 //origin[0] = 0.;
306 origin[1] = AliACORDEConstants::Instance()->Depth();
307 //origin[1] = 900;
308 origin[2] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
309 TMath::Cos(azimuth*kDegrad);
310 //origin[2] = 0.;
311
312 for (Int_t ix = 0; ix < fNx; ix++ ) {
313 for (Int_t iz = 0; iz < fNz; iz++ ) {
314 Rndm(random,6);
315 origin[0]+=ix*dx+2*(random[1]-0.5)*fOsigma[0];
316 origin[2]+=iz*dz+2*(random[2]-0.5)*fOsigma[2];
317 if ( random[4] < 0.5 ) {
318 origin[0] = -1*origin[0];
319 }
320 if ( random[5] < 0.5 ) {
321 origin[2] = -1*origin[2];
322 }
323
324 pmom = -(fPMin + random[3] *(fPMax - fPMax) ); // Always downwards
325 p[0] = TMath::Sin(zenith*kDegrad)*TMath::Sin(azimuth*kDegrad)*pmom;
326 p[1] = TMath::Cos(zenith*kDegrad)*pmom;
327 p[2] = TMath::Sin(zenith*kDegrad)*TMath::Cos(azimuth*kDegrad)*pmom;
328
329 PushTrack(fTrackIt, -1, fIpart, p, origin, polar, 0, kPPrimary, nt);
330 }
331
332 }
333
334}
335
336//____________________________________________________________________________
337void AliGenACORDE::SetGridRange(Int_t nx,Float_t xwidth, Int_t nz, Float_t zwidth)
338{
339 //
340 // Define the grid
341 // This data shuold be used for Muon bundles generation.
342 //
343 fXwidth=xwidth;
344 fNx=nx;
345 fZwidth=zwidth;
346 fNz=nz;
347
348 // Print a message about the use, if the Mode has not been set, or
349 // it has to a different Mode.
350 if ( fCRMode != kMuonBundle ) {
351 Warning("SetRange","You have been specified a grid to generate muon bundles, but seems that you haven't choose this generation mode, or you have already select a different one");
352 fMuonGrid = kTRUE;
353 }
354}
355
356//____________________________________________________________________________
357void AliGenACORDE::InitApWeightFactors()
358{
359 //
360 // This factors will be to correct the zenithal angle distribution
361 // acording the momentum
362
363 //
364 // Fill the array for the flux zenith angle dependence.
365 // at the index 0 of fAp[] will be the "factor" if we have a muon
366 // of 0 GeV.
367 Float_t a[6] = {-1.61, -1.50, -1.28, -0.94, -0.61, -0.22};
368 Float_t p[6] = { 0., 10., 30., 100., 300., 1000.};
369
370 // Get the information from the momentum
371 Int_t pEnd = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
372 // Initialize the Array of floats to hold the a(p) factors.
373 fAp = new TArrayF(pEnd);
374
375 Int_t index = 0;
376
377 for (Int_t i = 0; i < pEnd; i++ ) {
378 Float_t currentP = ((Float_t)i)*fPResolution;
379 if ( currentP < p[1] ) index = 0;
380 else if ( currentP >= p[1] && currentP < p[2] ) index = 1;
381 else if ( currentP >= p[2] && currentP < p[3] ) index = 2;
382 else if ( currentP >= p[3] && currentP < p[4] ) index = 3;
383 else if ( currentP >= p[4] ) index = 4;
384
385 Float_t ap = (currentP -p[index])*(a[index+1] - a[index])/
386 (p[index+1] - p[index]) + a[index];
387 fAp->AddAt(ap, i);
388 }
389
390}
391
392//___________________________________________________________________________
393void AliGenACORDE::InitMomentumGeneration()
394{
395 //
396 // Initialize a funtion to generate the momentum randomly
397 // acording this function.
398 //
399
400 // Check if we nned to initialize the function
401 if ( fPMin != fPMax ) {
402
403 // Check also if the function have been defined yet.
404 if ( !fMomentumDist ) {
405
406 // If not, use this function
407 const char* y = "log10(x)";
408
409 const char* h1Coef = "[0]*( %s*%s*%s/2 - (5*%s*%s/2) + 3*%s )";
410 const char* h2Coef = "[1]*( (-2*%s*%s*%s/3) + (3*%s*%s) - 10*%s/3 + 1 )";
411 const char* h3Coef = "[2]*( %s*%s*%s/6 - %s*%s/2 + %s/3 )";
412 const char* s2Coef = "[3]*( %s*%s*%s/3 - 2*%s*%s + 11*%s/3 - 2 )";
413
414 const char* h = "%s + %s + %s + %s";
415 const char* flux = "pow(10., %s)";
416 const char* normalizedFlux = "0.86*x*x*x*pow(10., %s)";
417 const char* paramNames[4] = {"H1", "H2", "H3", "S1"};
418
419 char buffer1[1024];
420 char buffer2[1024];
421 char buffer3[1024];
422 char buffer4[1024];
423 char buffer5[1024];
424 char buffer6[1024];
425 char buffer7[1024];
426
7e0f247e 427 snprintf(buffer1, 1023, h1Coef, y, y, y, y, y, y);
428 snprintf(buffer2, 1023, h2Coef, y, y, y, y, y, y);
429 snprintf(buffer3, 1023, h3Coef, y, y, y, y, y, y);
430 snprintf(buffer4, 1023, s2Coef, y, y, y, y, y, y);
b86e74f5 431
7e0f247e 432 snprintf(buffer5, 1023, h, buffer1, buffer2, buffer3, buffer4);
b86e74f5 433
7e0f247e 434 snprintf(buffer6, 1023, flux, buffer5);
b86e74f5 435
436 fMomentumDist = new TF1("fMomentumDist", buffer6, fPMin, fPMax);
7e0f247e 437 snprintf(buffer7, 1023, normalizedFlux, buffer5);
b86e74f5 438 fUnfoldedMomentumDist = new TF1("fUnfoldedMomentumDist", buffer7, fPMin, fPMax);
439 for (Int_t i = 0; i < 4; i++ ) {
440 fMomentumDist->SetParName(i, paramNames[i]);
441 fUnfoldedMomentumDist->SetParName(i, paramNames[i]);
442 }
443
444 fMomentumDist->SetParameter(0, 0.133);
445 fMomentumDist->SetParameter(1, -2.521);
446 fMomentumDist->SetParameter(2, -5.78);
447 fMomentumDist->SetParameter(3, -2.11);
448
449 fUnfoldedMomentumDist->SetParameter(0, 0.133);
450 fUnfoldedMomentumDist->SetParameter(1, -2.521);
451 fUnfoldedMomentumDist->SetParameter(2, -5.78);
452 fUnfoldedMomentumDist->SetParameter(3, -2.11);
453
454 }
455
456 }
457
458}
459
460//____________________________________________________________________________
461void AliGenACORDE::InitZenithalAngleGeneration()
462{
463 //
464 // Initalize a distribution function for the zenith angle.
465 // This angle will be obtained randomly acording this function.
466 // The generated angles will been in degrees.
467
468 // Check if we need to create the function.
469 if ( fZenithMin != fZenithMax ) {
470
471 // Check also if another function have been defined.
472 if ( !fZenithDist ) {
473
474 // initialize the momentum dependent coefficients, a(p)
475 this->InitApWeightFactors();
476
477 Int_t pEnd = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
478 char name[26];
479 char title[52];
480 fPDist = new TClonesArray("TH1F", pEnd);
481 TClonesArray &mom = *fPDist;
482 TH1F* zenith = 0;
483 Float_t weight = 0;
484 for ( Int_t i = 0; i < pEnd; i++ ) {
485 // Fill the distribution
7e0f247e 486 snprintf(name, 25, "zenith%d", i+1);
487 snprintf(title, 51, "Zenith distribution, p=%f", fPMin+(Float_t)i);
b86e74f5 488 zenith = new(mom[i]) TH1F(name, title, TMath::Abs(TMath::Nint(fZenithMax-fZenithMin)), TMath::Cos(fZenithMax*TMath::Pi()/180), TMath::Cos(fZenithMin*TMath::Pi()/180));
489
490 // Make a loop for the angle and fill the histogram for the weight
491 Int_t steps = 1000;
492 Float_t value = 0;
493 for (Int_t j = 0; j < steps; j++ ) {
494 value = TMath::Cos(fZenithMin*TMath::Pi()/180) + (Float_t)j * ( TMath::Cos(fZenithMax*TMath::Pi()/180) - TMath::Cos(fZenithMin*TMath::Pi()/180))/1000;
495 weight = 1 + fAp->At(i)*(1 - value);
496 zenith->Fill(value, weight);
497 }
498
499 }
500
501 }
502
503 }
504
505}
506
507//____________________________________________________________________________
508Float_t AliGenACORDE::GetZenithAngle(Float_t mom) const
509{
510
511 Float_t zenith = 0.;
512 // Check if you need to generate a constant zenith angle.
513 if ( !fZenithDist ) {
514 // Check if you have defined an array of momentum functions
515 if ( fPDist ) {
516 Int_t pIndex = TMath::Abs(TMath::Nint(mom));
517 TH1F* cosZenithAngle = (TH1F*)fPDist->UncheckedAt(pIndex);
518 Float_t tmpzenith = TMath::ACos(cosZenithAngle->GetRandom());
519 // Correct the value
520 zenith = kRaddeg*tmpzenith;
521 return zenith;
522 } else {
523
524 if ( fCRMode != kMuonFlux ) {
525 // If you aren't generating muons obeying any ditribution
526 // only generate a flat zenith angle, acording the input settings
527 Float_t random[2];
528 Rndm(random, 2);
529 zenith = fZenithMin + (fZenithMax - fZenithMin)*random[0];
530
531 } else {
532 // Even if you are generating muons acording some distribution,
533 // but you don't want to ...
534 zenith = fZenithMin;
535
536 }
537
538 }
539 } else {
540 zenith = fZenithDist->GetRandom();
541 }
542
543 return zenith;
544}
545
546//_____________________________________________________________________________
547Float_t AliGenACORDE::GetMomentum() const
548{
549 //
550 //
551 //
552 return fMomentumDist->GetRandom();
553}