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