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