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),
80 //_____________________________________________________________________________
81 AliGenACORDE::AliGenACORDE(Int_t npart)
82 : AliGenerator(npart),
84 fCRMode(kSingleMuons),
99 fUnfoldedMomentumDist(0),
108 fTitle = "Cosmic Muons generator";
110 // Set the origin above the vertex, on the surface.
112 fOrigin[1] = AliACORDEConstants::Instance()->Depth(); // At the surface by default.
116 //_____________________________________________________________________________
117 AliGenACORDE::~AliGenACORDE()
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;
129 //_____________________________________________________________________________
130 void AliGenACORDE::Generate()
133 // Generate on one trigger
134 // Call the respective method inside the loop for the number
135 // of tracks per trigger.
137 for (Int_t i = 0; i < fNParticles; i++ ) {
139 if ( fCRMode == kMuonBundle ) {
140 this->GenerateOneMuonBundle();
142 } else if ( fCRMode == kSingleMuons ) {
143 this->GenerateOneSingleMuon(kTRUE);
146 // Generate only single muons following the parametrizations
147 // for atmospheric muons.
148 this->GenerateOneSingleMuon();
155 //_____________________________________________________________________________
156 void AliGenACORDE::Init()
159 // Initialize some internal methods.
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");
168 // Determine some specific data members.
169 fPRange = TMath::Abs(fPMax-fPMin);
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!");
177 if((!TestBit(kPtRange))&&(!TestBit(kMomentumRange)))
178 Fatal("Init","You should set either the momentum or the pt range!");
180 } else if ( fCRMode == kMuonBundle ) {
181 fCRModeName = new TString("Muon Bundles");
183 } else if ( fCRMode == kMuonFlux ) {
184 fCRModeName = new TString("Muon Fluxes");
185 // Initialize the ditribution functions.
186 this->InitMomentumGeneration();
187 this->InitZenithalAngleGeneration();
190 Fatal("Generate", "Generation Mode unknown!\n");
196 //____________________________________________________________________________
197 void AliGenACORDE::GenerateOneSingleMuon(Bool_t withFlatMomentum)
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
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.
212 Float_t polar[3]= {0,0,0}; // Polarization parameters
219 // Take the azimuth random.
221 Float_t azimuth = fAzimuthMin + (fAzimuthMax-fAzimuthMin)*random[0];
222 Float_t zenith = fZenithMin + (fZenithMax - fZenithMin)*random[1];
224 if ( withFlatMomentum ) {
226 if(TestBit(kMomentumRange)) {
227 pmom = -( fPMin + random[0]*(fPMax - fPMin) ); // always downwards.
228 pt = pmom*TMath::Sin(zenith*kDegrad);
230 pt = -( fPtMin + random[1]*(fPtMax - fPtMin)); // always downwards.
231 pmom = pt/TMath::Sin(zenith*kDegrad);
235 if ( fMomentumDist ) {
236 pmom = -this->GetMomentum(); // Always downwards.
240 zenith = this->GetZenithAngle(pmom); // In degrees
241 pt = pmom*TMath::Sin(zenith*kDegrad);
244 p[0] = pt*TMath::Sin(azimuth*kDegrad);
245 p[1] = pmom*TMath::Cos(zenith*kDegrad);
246 p[2] = pt*TMath::Cos(azimuth*kDegrad);
248 // Finaly the origin, with the smearing
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]));
254 origin[1] = AliACORDEConstants::Instance()->Depth();
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]));
260 // Put the track on the stack.
261 PushTrack(fTrackIt,-1,fIpart,p,origin,polar,0,kPPrimary,nt);
265 //____________________________________________________________________________
266 void AliGenACORDE::GenerateOneMuonBundle()
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.
274 Float_t polar[3]= {0,0,0}; // Polarization parameters
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;
287 // Generate the kinematics a la AliGenScan (Andreas Morchs)
303 origin[0] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
304 TMath::Sin(azimuth*kDegrad);
306 origin[1] = AliACORDEConstants::Instance()->Depth();
308 origin[2] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
309 TMath::Cos(azimuth*kDegrad);
312 for (Int_t ix = 0; ix < fNx; ix++ ) {
313 for (Int_t iz = 0; iz < fNz; iz++ ) {
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];
320 if ( random[5] < 0.5 ) {
321 origin[2] = -1*origin[2];
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;
329 PushTrack(fTrackIt, -1, fIpart, p, origin, polar, 0, kPPrimary, nt);
336 //____________________________________________________________________________
337 void AliGenACORDE::SetGridRange(Int_t nx,Float_t xwidth, Int_t nz, Float_t zwidth)
341 // This data shuold be used for Muon bundles generation.
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");
356 //____________________________________________________________________________
357 void AliGenACORDE::InitApWeightFactors()
360 // This factors will be to correct the zenithal angle distribution
361 // acording the momentum
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
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.};
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);
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;
385 Float_t ap = (currentP -p[index])*(a[index+1] - a[index])/
386 (p[index+1] - p[index]) + a[index];
392 //___________________________________________________________________________
393 void AliGenACORDE::InitMomentumGeneration()
396 // Initialize a funtion to generate the momentum randomly
397 // acording this function.
400 // Check if we nned to initialize the function
401 if ( fPMin != fPMax ) {
403 // Check also if the function have been defined yet.
404 if ( !fMomentumDist ) {
406 // If not, use this function
407 const char* y = "log10(x)";
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 )";
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"};
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);
432 snprintf(buffer5, 1023, h, buffer1, buffer2, buffer3, buffer4);
434 snprintf(buffer6, 1023, flux, buffer5);
436 fMomentumDist = new TF1("fMomentumDist", buffer6, fPMin, fPMax);
437 snprintf(buffer7, 1023, normalizedFlux, buffer5);
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]);
444 fMomentumDist->SetParameter(0, 0.133);
445 fMomentumDist->SetParameter(1, -2.521);
446 fMomentumDist->SetParameter(2, -5.78);
447 fMomentumDist->SetParameter(3, -2.11);
449 fUnfoldedMomentumDist->SetParameter(0, 0.133);
450 fUnfoldedMomentumDist->SetParameter(1, -2.521);
451 fUnfoldedMomentumDist->SetParameter(2, -5.78);
452 fUnfoldedMomentumDist->SetParameter(3, -2.11);
460 //____________________________________________________________________________
461 void AliGenACORDE::InitZenithalAngleGeneration()
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.
468 // Check if we need to create the function.
469 if ( fZenithMin != fZenithMax ) {
471 // Check also if another function have been defined.
472 if ( !fZenithDist ) {
474 // initialize the momentum dependent coefficients, a(p)
475 this->InitApWeightFactors();
477 Int_t pEnd = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
480 fPDist = new TClonesArray("TH1F", pEnd);
481 TClonesArray &mom = *fPDist;
484 for ( Int_t i = 0; i < pEnd; i++ ) {
485 // Fill the distribution
486 snprintf(name, 25, "zenith%d", i+1);
487 snprintf(title, 51, "Zenith distribution, p=%f", fPMin+(Float_t)i);
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));
490 // Make a loop for the angle and fill the histogram for the weight
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);
507 //____________________________________________________________________________
508 Float_t AliGenACORDE::GetZenithAngle(Float_t mom) const
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
516 Int_t pIndex = TMath::Abs(TMath::Nint(mom));
517 TH1F* cosZenithAngle = (TH1F*)fPDist->UncheckedAt(pIndex);
518 Float_t tmpzenith = TMath::ACos(cosZenithAngle->GetRandom());
520 zenith = kRaddeg*tmpzenith;
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
529 zenith = fZenithMin + (fZenithMax - fZenithMin)*random[0];
532 // Even if you are generating muons acording some distribution,
533 // but you don't want to ...
540 zenith = fZenithDist->GetRandom();
546 //_____________________________________________________________________________
547 Float_t AliGenACORDE::GetMomentum() const
552 return fMomentumDist->GetRandom();