Line data Source code
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 :
49 12 : ClassImp(AliGenACORDE)
50 :
51 : //_____________________________________________________________________________
52 : AliGenACORDE::AliGenACORDE()
53 0 : : AliGenerator(),
54 0 : fIpart(0),
55 0 : fCRMode(kSingleMuons),
56 0 : fCRModeName(0),
57 0 : fXwidth(0),
58 0 : fNx(1),
59 0 : fZwidth(0),
60 0 : fNz(1),
61 0 : fMuonGrid(kFALSE),
62 0 : fZenithMin(0),
63 0 : fZenithMax(0),
64 0 : fAzimuthMin(0),
65 0 : fAzimuthMax(0),
66 0 : fPRange(0),
67 0 : fPResolution(1),
68 0 : fAp(0),
69 0 : fMomentumDist(0),
70 0 : fUnfoldedMomentumDist(0),
71 0 : fZenithDist(0),
72 0 : fPDist(0),
73 0 : fNParticles(0)
74 0 : {
75 : //
76 : // Default ctor.
77 : //
78 0 : }
79 :
80 : //_____________________________________________________________________________
81 : AliGenACORDE::AliGenACORDE(Int_t npart)
82 0 : : AliGenerator(npart),
83 0 : fIpart(kMuonMinus),
84 0 : fCRMode(kSingleMuons),
85 0 : fCRModeName(0),
86 0 : fXwidth(0),
87 0 : fNx(1),
88 0 : fZwidth(0),
89 0 : fNz(1),
90 0 : fMuonGrid(kFALSE),
91 0 : fZenithMin(0),
92 0 : fZenithMax(0),
93 0 : fAzimuthMin(0),
94 0 : fAzimuthMax(0),
95 0 : fPRange(0),
96 0 : fPResolution(1),
97 0 : fAp(0),
98 0 : fMomentumDist(0),
99 0 : fUnfoldedMomentumDist(0),
100 0 : fZenithDist(0),
101 0 : fPDist(0),
102 0 : fNParticles(0)
103 0 : {
104 : //
105 : // Standard ctor.
106 : //
107 0 : fName = "ACORDE";
108 0 : fTitle = "Cosmic Muons generator";
109 :
110 : // Set the origin above the vertex, on the surface.
111 0 : fOrigin[0] = 0.;
112 0 : fOrigin[1] = AliACORDEConstants::Instance()->Depth(); // At the surface by default.
113 0 : fOrigin[2] = 0.;
114 0 : }
115 :
116 : //_____________________________________________________________________________
117 0 : AliGenACORDE::~AliGenACORDE()
118 0 : {
119 : //
120 : // Default dtor.
121 : //
122 0 : if ( fPDist ) {fPDist->Delete(); delete fPDist; fPDist = 0;}
123 0 : if ( fUnfoldedMomentumDist ) delete fUnfoldedMomentumDist;
124 0 : if ( fMomentumDist ) delete fMomentumDist;
125 0 : if ( fAp ) delete fAp;
126 0 : if ( fCRModeName ) delete fCRModeName;
127 0 : }
128 :
129 : //_____________________________________________________________________________
130 : void 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 :
137 0 : for (Int_t i = 0; i < fNParticles; i++ ) {
138 :
139 0 : if ( fCRMode == kMuonBundle ) {
140 0 : this->GenerateOneMuonBundle();
141 :
142 0 : } else if ( fCRMode == kSingleMuons ) {
143 0 : this->GenerateOneSingleMuon(kTRUE);
144 :
145 0 : } else {
146 : // Generate only single muons following the parametrizations
147 : // for atmospheric muons.
148 0 : this->GenerateOneSingleMuon();
149 :
150 : }
151 :
152 : }
153 0 : }
154 :
155 : //_____________________________________________________________________________
156 : void AliGenACORDE::Init()
157 : {
158 : //
159 : // Initialize some internal methods.
160 : //
161 :
162 :
163 0 : printf("**************************************************************\n");
164 0 : printf("<<< *** Starting the AliGenACORDE cosmic generator ******** >>>\n");
165 0 : printf("<<< *** No. of muons generated at the surface of P2: %d, * >>>\n",fNParticles);
166 0 : printf("**************************************************************\n");
167 :
168 : // Determine some specific data members.
169 0 : fPRange = TMath::Abs(fPMax-fPMin);
170 :
171 0 : if ( fCRMode == kSingleMuons ) {
172 0 : fCRModeName = new TString("Single Muons");
173 : // Initialisation, check consistency of selected ranges
174 0 : if(TestBit(kPtRange)&&TestBit(kMomentumRange))
175 0 : Fatal("Init","You should not set the momentum range and the pt range!");
176 :
177 0 : if((!TestBit(kPtRange))&&(!TestBit(kMomentumRange)))
178 0 : Fatal("Init","You should set either the momentum or the pt range!");
179 :
180 0 : } else if ( fCRMode == kMuonBundle ) {
181 0 : fCRModeName = new TString("Muon Bundles");
182 :
183 0 : } else if ( fCRMode == kMuonFlux ) {
184 0 : fCRModeName = new TString("Muon Fluxes");
185 : // Initialize the ditribution functions.
186 0 : this->InitMomentumGeneration();
187 0 : this->InitZenithalAngleGeneration();
188 :
189 0 : } else {
190 0 : Fatal("Generate", "Generation Mode unknown!\n");
191 :
192 : }
193 :
194 0 : }
195 :
196 : //____________________________________________________________________________
197 : void 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 0 : Float_t polar[3]= {0,0,0}; // Polarization parameters
213 0 : Float_t origin[3];
214 0 : Int_t nt;
215 0 : Float_t p[3];
216 : Float_t pmom, pt;
217 0 : Float_t random[6];
218 :
219 : // Take the azimuth random.
220 0 : Rndm(random, 2);
221 0 : Float_t azimuth = fAzimuthMin + (fAzimuthMax-fAzimuthMin)*random[0];
222 0 : Float_t zenith = fZenithMin + (fZenithMax - fZenithMin)*random[1];
223 :
224 0 : if ( withFlatMomentum ) {
225 0 : Rndm(random,3);
226 0 : if(TestBit(kMomentumRange)) {
227 0 : pmom = -( fPMin + random[0]*(fPMax - fPMin) ); // always downwards.
228 0 : pt = pmom*TMath::Sin(zenith*kDegrad);
229 0 : } else {
230 0 : pt = -( fPtMin + random[1]*(fPtMax - fPtMin)); // always downwards.
231 0 : pmom = pt/TMath::Sin(zenith*kDegrad);
232 : }
233 :
234 : } else {
235 0 : if ( fMomentumDist ) {
236 0 : pmom = -this->GetMomentum(); // Always downwards.
237 0 : } else {
238 0 : pmom = -fPMin;
239 : }
240 0 : zenith = this->GetZenithAngle(pmom); // In degrees
241 0 : pt = pmom*TMath::Sin(zenith*kDegrad);
242 : }
243 :
244 0 : p[0] = pt*TMath::Sin(azimuth*kDegrad);
245 0 : p[1] = pmom*TMath::Cos(zenith*kDegrad);
246 0 : p[2] = pt*TMath::Cos(azimuth*kDegrad);
247 :
248 : // Finaly the origin, with the smearing
249 0 : Rndm(random,6);
250 0 : origin[0] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
251 0 : TMath::Sin(azimuth*kDegrad)
252 0 : + fOsigma[0]* TMath::Cos(2*random[0]*TMath::Pi())*TMath::Sqrt(-2*TMath::Log(random[1]));
253 :
254 0 : origin[1] = AliACORDEConstants::Instance()->Depth();
255 :
256 0 : origin[2] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
257 0 : TMath::Cos(azimuth*kDegrad)
258 0 : + 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 0 : PushTrack(fTrackIt,-1,fIpart,p,origin,polar,0,kPPrimary,nt);
262 :
263 0 : }
264 :
265 : //____________________________________________________________________________
266 : void 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 0 : Float_t polar[3]= {0,0,0}; // Polarization parameters
275 0 : Float_t origin[3];
276 0 : Float_t p[3];
277 0 : Int_t nt;
278 : Float_t pmom;
279 0 : Float_t random[6];
280 :
281 0 : Rndm(random, 3);
282 0 : Float_t zenith = fZenithMin + (fZenithMax - fZenithMin)*random[1];
283 0 : 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 0 : if ( fNx > 0 ) {
290 0 : dx = fXwidth/fNx;
291 0 : } else {
292 : dx = 1e10;
293 : //dx = 100.;
294 : }
295 :
296 0 : if ( fNz > 0 ) {
297 0 : dz = fZwidth/fNz;
298 0 : } else {
299 : dz = 1e10;
300 : //dz = 100.;
301 : }
302 :
303 0 : origin[0] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
304 0 : TMath::Sin(azimuth*kDegrad);
305 : //origin[0] = 0.;
306 0 : origin[1] = AliACORDEConstants::Instance()->Depth();
307 : //origin[1] = 900;
308 0 : origin[2] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
309 0 : TMath::Cos(azimuth*kDegrad);
310 : //origin[2] = 0.;
311 :
312 0 : for (Int_t ix = 0; ix < fNx; ix++ ) {
313 0 : for (Int_t iz = 0; iz < fNz; iz++ ) {
314 0 : Rndm(random,6);
315 0 : origin[0]+=ix*dx+2*(random[1]-0.5)*fOsigma[0];
316 0 : origin[2]+=iz*dz+2*(random[2]-0.5)*fOsigma[2];
317 0 : if ( random[4] < 0.5 ) {
318 0 : origin[0] = -1*origin[0];
319 0 : }
320 0 : if ( random[5] < 0.5 ) {
321 0 : origin[2] = -1*origin[2];
322 0 : }
323 :
324 0 : pmom = -(fPMin + random[3] *(fPMax - fPMax) ); // Always downwards
325 0 : p[0] = TMath::Sin(zenith*kDegrad)*TMath::Sin(azimuth*kDegrad)*pmom;
326 0 : p[1] = TMath::Cos(zenith*kDegrad)*pmom;
327 0 : p[2] = TMath::Sin(zenith*kDegrad)*TMath::Cos(azimuth*kDegrad)*pmom;
328 :
329 0 : PushTrack(fTrackIt, -1, fIpart, p, origin, polar, 0, kPPrimary, nt);
330 : }
331 :
332 : }
333 :
334 0 : }
335 :
336 : //____________________________________________________________________________
337 : void 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 0 : fXwidth=xwidth;
344 0 : fNx=nx;
345 0 : fZwidth=zwidth;
346 0 : 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 0 : if ( fCRMode != kMuonBundle ) {
351 0 : 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 0 : fMuonGrid = kTRUE;
353 0 : }
354 0 : }
355 :
356 : //____________________________________________________________________________
357 : void 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 0 : Float_t a[6] = {-1.61, -1.50, -1.28, -0.94, -0.61, -0.22};
368 0 : Float_t p[6] = { 0., 10., 30., 100., 300., 1000.};
369 :
370 : // Get the information from the momentum
371 0 : Int_t pEnd = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
372 : // Initialize the Array of floats to hold the a(p) factors.
373 0 : fAp = new TArrayF(pEnd);
374 :
375 : Int_t index = 0;
376 :
377 0 : for (Int_t i = 0; i < pEnd; i++ ) {
378 0 : Float_t currentP = ((Float_t)i)*fPResolution;
379 0 : if ( currentP < p[1] ) index = 0;
380 0 : else if ( currentP >= p[1] && currentP < p[2] ) index = 1;
381 0 : else if ( currentP >= p[2] && currentP < p[3] ) index = 2;
382 0 : else if ( currentP >= p[3] && currentP < p[4] ) index = 3;
383 0 : else if ( currentP >= p[4] ) index = 4;
384 :
385 0 : Float_t ap = (currentP -p[index])*(a[index+1] - a[index])/
386 0 : (p[index+1] - p[index]) + a[index];
387 0 : fAp->AddAt(ap, i);
388 : }
389 :
390 0 : }
391 :
392 : //___________________________________________________________________________
393 : void 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 0 : if ( fPMin != fPMax ) {
402 :
403 : // Check also if the function have been defined yet.
404 0 : 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 0 : const char* paramNames[4] = {"H1", "H2", "H3", "S1"};
418 :
419 0 : char buffer1[1024];
420 0 : char buffer2[1024];
421 0 : char buffer3[1024];
422 0 : char buffer4[1024];
423 0 : char buffer5[1024];
424 0 : char buffer6[1024];
425 0 : char buffer7[1024];
426 :
427 0 : snprintf(buffer1, 1023, h1Coef, y, y, y, y, y, y);
428 0 : snprintf(buffer2, 1023, h2Coef, y, y, y, y, y, y);
429 0 : snprintf(buffer3, 1023, h3Coef, y, y, y, y, y, y);
430 0 : snprintf(buffer4, 1023, s2Coef, y, y, y, y, y, y);
431 :
432 0 : snprintf(buffer5, 1023, h, buffer1, buffer2, buffer3, buffer4);
433 :
434 0 : snprintf(buffer6, 1023, flux, buffer5);
435 :
436 0 : fMomentumDist = new TF1("fMomentumDist", buffer6, fPMin, fPMax);
437 0 : snprintf(buffer7, 1023, normalizedFlux, buffer5);
438 0 : fUnfoldedMomentumDist = new TF1("fUnfoldedMomentumDist", buffer7, fPMin, fPMax);
439 0 : for (Int_t i = 0; i < 4; i++ ) {
440 0 : fMomentumDist->SetParName(i, paramNames[i]);
441 0 : fUnfoldedMomentumDist->SetParName(i, paramNames[i]);
442 : }
443 :
444 0 : fMomentumDist->SetParameter(0, 0.133);
445 0 : fMomentumDist->SetParameter(1, -2.521);
446 0 : fMomentumDist->SetParameter(2, -5.78);
447 0 : fMomentumDist->SetParameter(3, -2.11);
448 :
449 0 : fUnfoldedMomentumDist->SetParameter(0, 0.133);
450 0 : fUnfoldedMomentumDist->SetParameter(1, -2.521);
451 0 : fUnfoldedMomentumDist->SetParameter(2, -5.78);
452 0 : fUnfoldedMomentumDist->SetParameter(3, -2.11);
453 :
454 0 : }
455 :
456 : }
457 :
458 0 : }
459 :
460 : //____________________________________________________________________________
461 : void 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 0 : if ( fZenithMin != fZenithMax ) {
470 :
471 : // Check also if another function have been defined.
472 0 : if ( !fZenithDist ) {
473 :
474 : // initialize the momentum dependent coefficients, a(p)
475 0 : this->InitApWeightFactors();
476 :
477 0 : Int_t pEnd = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
478 0 : char name[26];
479 0 : char title[52];
480 0 : fPDist = new TClonesArray("TH1F", pEnd);
481 : TClonesArray &mom = *fPDist;
482 : TH1F* zenith = 0;
483 : Float_t weight = 0;
484 0 : for ( Int_t i = 0; i < pEnd; i++ ) {
485 : // Fill the distribution
486 0 : snprintf(name, 25, "zenith%d", i+1);
487 0 : snprintf(title, 51, "Zenith distribution, p=%f", fPMin+(Float_t)i);
488 0 : 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 0 : for (Int_t j = 0; j < steps; j++ ) {
494 0 : value = TMath::Cos(fZenithMin*TMath::Pi()/180) + (Float_t)j * ( TMath::Cos(fZenithMax*TMath::Pi()/180) - TMath::Cos(fZenithMin*TMath::Pi()/180))/1000;
495 0 : weight = 1 + fAp->At(i)*(1 - value);
496 0 : zenith->Fill(value, weight);
497 : }
498 :
499 : }
500 :
501 0 : }
502 :
503 : }
504 :
505 0 : }
506 :
507 : //____________________________________________________________________________
508 : Float_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 0 : if ( !fZenithDist ) {
514 : // Check if you have defined an array of momentum functions
515 0 : if ( fPDist ) {
516 0 : Int_t pIndex = TMath::Abs(TMath::Nint(mom));
517 0 : TH1F* cosZenithAngle = (TH1F*)fPDist->UncheckedAt(pIndex);
518 0 : Float_t tmpzenith = TMath::ACos(cosZenithAngle->GetRandom());
519 : // Correct the value
520 0 : zenith = kRaddeg*tmpzenith;
521 : return zenith;
522 : } else {
523 :
524 0 : 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 0 : Float_t random[2];
528 0 : Rndm(random, 2);
529 0 : zenith = fZenithMin + (fZenithMax - fZenithMin)*random[0];
530 :
531 0 : } else {
532 : // Even if you are generating muons acording some distribution,
533 : // but you don't want to ...
534 0 : zenith = fZenithMin;
535 :
536 : }
537 :
538 : }
539 : } else {
540 0 : zenith = fZenithDist->GetRandom();
541 : }
542 :
543 0 : return zenith;
544 0 : }
545 :
546 : //_____________________________________________________________________________
547 : Float_t AliGenACORDE::GetMomentum() const
548 : {
549 : //
550 : //
551 : //
552 0 : return fMomentumDist->GetRandom();
553 : }
|
