1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 /////////////////////////////////////////////////////////////////////////////
20 // Contain parametrizations to generate atmospheric muons, and also
21 // to generate single muons and muon bundles at surface level.
25 <img src="picts/AliGenACORDEClass.gif">
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>.
36 /////////////////////////////////////////////////////////////////////////////
38 #include "AliGenACORDE.h"
40 #include <TMCProcess.h>
42 #include <TClonesArray.h>
49 ClassImp(AliGenACORDE)
51 //_____________________________________________________________________________
52 AliGenACORDE::AliGenACORDE()
55 fCRMode(kSingleMuons),
70 fUnfoldedMomentumDist(0),
79 //_____________________________________________________________________________
80 AliGenACORDE::AliGenACORDE(Int_t npart)
81 : AliGenerator(npart),
83 fCRMode(kSingleMuons),
98 fUnfoldedMomentumDist(0),
106 fTitle = "Cosmic Muons generator";
108 // Set the origin above the vertex, on the surface.
110 fOrigin[1] = AliACORDEConstants::Instance()->Depth(); // At the surface by default.
114 //_____________________________________________________________________________
115 AliGenACORDE::~AliGenACORDE()
120 if ( fPDist ) {fPDist->Delete(); delete fPDist; fPDist = 0;}
121 if ( fUnfoldedMomentumDist ) delete fUnfoldedMomentumDist;
122 if ( fMomentumDist ) delete fMomentumDist;
123 if ( fAp ) delete fAp;
124 if ( fCRModeName ) delete fCRModeName;
127 //_____________________________________________________________________________
128 void AliGenACORDE::Generate()
131 // Generate on one trigger
132 // Call the respective method inside the loop for the number
133 // of tracks per trigger.
135 for (Int_t i = 0; i < fNpart; i++ ) {
137 if ( fCRMode == kMuonBundle ) {
138 this->GenerateOneMuonBundle();
140 } else if ( fCRMode == kSingleMuons ) {
141 this->GenerateOneSingleMuon(kTRUE);
144 // Generate only single muons following the parametrizations
145 // for atmospheric muons.
146 this->GenerateOneSingleMuon();
153 //_____________________________________________________________________________
154 void AliGenACORDE::Init()
157 // Initialize some internal methods.
160 // Determine some specific data members.
161 fPRange = TMath::Abs(fPMax-fPMin);
163 if ( fCRMode == kSingleMuons ) {
164 fCRModeName = new TString("Single Muons");
165 // Initialisation, check consistency of selected ranges
166 if(TestBit(kPtRange)&&TestBit(kMomentumRange))
167 Fatal("Init","You should not set the momentum range and the pt range!");
169 if((!TestBit(kPtRange))&&(!TestBit(kMomentumRange)))
170 Fatal("Init","You should set either the momentum or the pt range!");
172 } else if ( fCRMode == kMuonBundle ) {
173 fCRModeName = new TString("Muon Bundles");
175 } else if ( fCRMode == kMuonFlux ) {
176 fCRModeName = new TString("Muon Fluxes");
177 // Initialize the ditribution functions.
178 this->InitMomentumGeneration();
179 this->InitZenithalAngleGeneration();
182 Fatal("Generate", "Generation Mode unknown!\n");
188 //____________________________________________________________________________
189 void AliGenACORDE::GenerateOneSingleMuon(Bool_t withFlatMomentum)
192 // Generate One Single Muon
193 // This method will generate only one muon.
194 // The momentum will be randomly flat distributed if
195 // the paremeter "withFlatMomentum" is set to kTRUE,
196 // otherwise the momentum will generate acordingly the parametrization
198 // and adpted from Bruno Alessandro's implementation with the
199 // CERNLIB to AliRoot.
200 // The "withFlatMomentum" parameter also will be used to generate
201 // the muons with a flat Zenithal angle distribution.
202 // Do the smearing here, so that means per track.
204 Float_t polar[3]= {0,0,0}; // Polarization parameters
211 // Take the azimuth random.
213 Float_t azimuth = fAzimuthMin + (fAzimuthMax-fAzimuthMin)*random[0];
214 Float_t zenith = fZenithMin + (fZenithMax - fZenithMin)*random[1];
216 if ( withFlatMomentum ) {
218 if(TestBit(kMomentumRange)) {
219 pmom = -( fPMin + random[0]*(fPMax - fPMin) ); // always downwards.
220 pt = pmom*TMath::Sin(zenith*kDegrad);
222 pt = -( fPtMin + random[1]*(fPtMax - fPtMin)); // always downwards.
223 pmom = pt/TMath::Sin(zenith*kDegrad);
227 if ( fMomentumDist ) {
228 pmom = -this->GetMomentum(); // Always downwards.
232 zenith = this->GetZenithAngle(pmom); // In degrees
233 pt = pmom*TMath::Sin(zenith*kDegrad);
236 p[0] = pt*TMath::Sin(azimuth*kDegrad);
237 p[1] = pmom*TMath::Cos(zenith*kDegrad);
238 p[2] = pt*TMath::Cos(azimuth*kDegrad);
240 // Finaly the origin, with the smearing
242 origin[0] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
243 TMath::Sin(azimuth*kDegrad)
244 + fOsigma[0]* TMath::Cos(2*random[0]*TMath::Pi())*TMath::Sqrt(-2*TMath::Log(random[1]));
246 origin[1] = AliACORDEConstants::Instance()->Depth();
248 origin[2] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
249 TMath::Cos(azimuth*kDegrad)
250 + fOsigma[2]* TMath::Cos(2*random[2]*TMath::Pi())*TMath::Sqrt(-2*TMath::Log(random[3]));
252 // Put the track on the stack.
253 PushTrack(fTrackIt,-1,fIpart,p,origin,polar,0,kPPrimary,nt);
257 //____________________________________________________________________________
258 void AliGenACORDE::GenerateOneMuonBundle()
261 // Generate One Muon Bundle method
262 // This method will generate a bunch of muons following the
263 // procedure of the AliGenScan class.
264 // These muons will be generated with flat momentum.
266 Float_t polar[3]= {0,0,0}; // Polarization parameters
274 Float_t zenith = fZenithMin + (fZenithMax - fZenithMin)*random[1];
275 Float_t azimuth = fAzimuthMin + (fAzimuthMax-fAzimuthMin)*random[2];
276 //Float_t zenith = 10;
277 //Float_t azimuth = 30;
279 // Generate the kinematics a la AliGenScan (Andreas Morchs)
295 origin[0] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
296 TMath::Sin(azimuth*kDegrad);
298 origin[1] = AliACORDEConstants::Instance()->Depth();
300 origin[2] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
301 TMath::Cos(azimuth*kDegrad);
304 for (Int_t ix = 0; ix < fNx; ix++ ) {
305 for (Int_t iz = 0; iz < fNz; iz++ ) {
307 origin[0]+=ix*dx+2*(random[1]-0.5)*fOsigma[0];
308 origin[2]+=iz*dz+2*(random[2]-0.5)*fOsigma[2];
309 if ( random[4] < 0.5 ) {
310 origin[0] = -1*origin[0];
312 if ( random[5] < 0.5 ) {
313 origin[2] = -1*origin[2];
316 pmom = -(fPMin + random[3] *(fPMax - fPMax) ); // Always downwards
317 p[0] = TMath::Sin(zenith*kDegrad)*TMath::Sin(azimuth*kDegrad)*pmom;
318 p[1] = TMath::Cos(zenith*kDegrad)*pmom;
319 p[2] = TMath::Sin(zenith*kDegrad)*TMath::Cos(azimuth*kDegrad)*pmom;
321 PushTrack(fTrackIt, -1, fIpart, p, origin, polar, 0, kPPrimary, nt);
328 //____________________________________________________________________________
329 void AliGenACORDE::SetGridRange(Int_t nx,Float_t xwidth, Int_t nz, Float_t zwidth)
333 // This data shuold be used for Muon bundles generation.
340 // Print a message about the use, if the Mode has not been set, or
341 // it has to a different Mode.
342 if ( fCRMode != kMuonBundle ) {
343 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");
348 //____________________________________________________________________________
349 void AliGenACORDE::InitApWeightFactors()
352 // This factors will be to correct the zenithal angle distribution
353 // acording the momentum
356 // Fill the array for the flux zenith angle dependence.
357 // at the index 0 of fAp[] will be the "factor" if we have a muon
359 Float_t a[6] = {-1.61, -1.50, -1.28, -0.94, -0.61, -0.22};
360 Float_t p[6] = { 0., 10., 30., 100., 300., 1000.};
362 // Get the information from the momentum
363 Int_t pEnd = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
364 // Initialize the Array of floats to hold the a(p) factors.
365 fAp = new TArrayF(pEnd);
369 for (Int_t i = 0; i < pEnd; i++ ) {
370 Float_t currentP = ((Float_t)i)*fPResolution;
371 if ( currentP < p[1] ) index = 0;
372 else if ( currentP >= p[1] && currentP < p[2] ) index = 1;
373 else if ( currentP >= p[2] && currentP < p[3] ) index = 2;
374 else if ( currentP >= p[3] && currentP < p[4] ) index = 3;
375 else if ( currentP >= p[4] ) index = 4;
377 Float_t ap = (currentP -p[index])*(a[index+1] - a[index])/
378 (p[index+1] - p[index]) + a[index];
384 //___________________________________________________________________________
385 void AliGenACORDE::InitMomentumGeneration()
388 // Initialize a funtion to generate the momentum randomly
389 // acording this function.
392 // Check if we nned to initialize the function
393 if ( fPMin != fPMax ) {
395 // Check also if the function have been defined yet.
396 if ( !fMomentumDist ) {
398 // If not, use this function
399 const char* y = "log10(x)";
401 const char* h1Coef = "[0]*( %s*%s*%s/2 - (5*%s*%s/2) + 3*%s )";
402 const char* h2Coef = "[1]*( (-2*%s*%s*%s/3) + (3*%s*%s) - 10*%s/3 + 1 )";
403 const char* h3Coef = "[2]*( %s*%s*%s/6 - %s*%s/2 + %s/3 )";
404 const char* s2Coef = "[3]*( %s*%s*%s/3 - 2*%s*%s + 11*%s/3 - 2 )";
406 const char* h = "%s + %s + %s + %s";
407 const char* flux = "pow(10., %s)";
408 const char* normalizedFlux = "0.86*x*x*x*pow(10., %s)";
409 const char* paramNames[4] = {"H1", "H2", "H3", "S1"};
419 sprintf(buffer1, h1Coef, y, y, y, y, y, y);
420 sprintf(buffer2, h2Coef, y, y, y, y, y, y);
421 sprintf(buffer3, h3Coef, y, y, y, y, y, y);
422 sprintf(buffer4, s2Coef, y, y, y, y, y, y);
424 sprintf(buffer5, h, buffer1, buffer2, buffer3, buffer4);
426 sprintf(buffer6, flux, buffer5);
428 fMomentumDist = new TF1("fMomentumDist", buffer6, fPMin, fPMax);
429 sprintf(buffer7, normalizedFlux, buffer5);
430 fUnfoldedMomentumDist = new TF1("fUnfoldedMomentumDist", buffer7, fPMin, fPMax);
431 for (Int_t i = 0; i < 4; i++ ) {
432 fMomentumDist->SetParName(i, paramNames[i]);
433 fUnfoldedMomentumDist->SetParName(i, paramNames[i]);
436 fMomentumDist->SetParameter(0, 0.133);
437 fMomentumDist->SetParameter(1, -2.521);
438 fMomentumDist->SetParameter(2, -5.78);
439 fMomentumDist->SetParameter(3, -2.11);
441 fUnfoldedMomentumDist->SetParameter(0, 0.133);
442 fUnfoldedMomentumDist->SetParameter(1, -2.521);
443 fUnfoldedMomentumDist->SetParameter(2, -5.78);
444 fUnfoldedMomentumDist->SetParameter(3, -2.11);
452 //____________________________________________________________________________
453 void AliGenACORDE::InitZenithalAngleGeneration()
456 // Initalize a distribution function for the zenith angle.
457 // This angle will be obtained randomly acording this function.
458 // The generated angles will been in degrees.
460 // Check if we need to create the function.
461 if ( fZenithMin != fZenithMax ) {
463 // Check also if another function have been defined.
464 if ( !fZenithDist ) {
466 // initialize the momentum dependent coefficients, a(p)
467 this->InitApWeightFactors();
469 Int_t pEnd = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
472 fPDist = new TClonesArray("TH1F", pEnd);
473 TClonesArray &mom = *fPDist;
476 for ( Int_t i = 0; i < pEnd; i++ ) {
477 // Fill the distribution
478 sprintf(name, "zenith%d", i+1);
479 sprintf(title, "Zenith distribution, p=%f", fPMin+(Float_t)i);
480 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));
482 // Make a loop for the angle and fill the histogram for the weight
485 for (Int_t j = 0; j < steps; j++ ) {
486 value = TMath::Cos(fZenithMin*TMath::Pi()/180) + (Float_t)j * ( TMath::Cos(fZenithMax*TMath::Pi()/180) - TMath::Cos(fZenithMin*TMath::Pi()/180))/1000;
487 weight = 1 + fAp->At(i)*(1 - value);
488 zenith->Fill(value, weight);
499 //____________________________________________________________________________
500 Float_t AliGenACORDE::GetZenithAngle(Float_t mom) const
504 // Check if you need to generate a constant zenith angle.
505 if ( !fZenithDist ) {
506 // Check if you have defined an array of momentum functions
508 Int_t pIndex = TMath::Abs(TMath::Nint(mom));
509 TH1F* cosZenithAngle = (TH1F*)fPDist->UncheckedAt(pIndex);
510 Float_t tmpzenith = TMath::ACos(cosZenithAngle->GetRandom());
512 zenith = kRaddeg*tmpzenith;
516 if ( fCRMode != kMuonFlux ) {
517 // If you aren't generating muons obeying any ditribution
518 // only generate a flat zenith angle, acording the input settings
521 zenith = fZenithMin + (fZenithMax - fZenithMin)*random[0];
524 // Even if you are generating muons acording some distribution,
525 // but you don't want to ...
532 zenith = fZenithDist->GetRandom();
538 //_____________________________________________________________________________
539 Float_t AliGenACORDE::GetMomentum() const
544 return fMomentumDist->GetRandom();