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 Revision 1.4 2002/11/21 19:34:52 alibrary
19 Removing AliMC and AliMCProcess
21 Revision 1.3 2002/10/23 06:47:56 alibrary
22 Introducing Riostream.h
24 Revision 1.2 2002/10/14 14:55:34 hristov
25 Merging the VirtualMC branch to the main development branch (HEAD)
27 Revision 1.1.2.1 2002/10/10 14:40:31 hristov
28 Updating VirtualMC to v3-09-02
30 Revision 1.1 2002/10/07 11:25:28 gamez
31 First version, generation of cosmic muons on the surface
36 /////////////////////////////////////////////////////////////////////////////
38 // Contain parametrizations to generate atmospheric muons, and also
39 // to generate single muons and muon bundles at surface level.
43 <img src="picts/AliGenCRTClass.gif">
46 <font size=+2 color=red>
47 <p>The responsible person for this module is
48 <a href="mailto:Enrique.Gamez.Flores@cern.ch">Enrique Gamez</a>.
54 /////////////////////////////////////////////////////////////////////////////
56 #include <Riostream.h>
57 #include <TMCProcess.h>
62 #include "AliGenCRT.h"
67 //_____________________________________________________________________________
68 AliGenCRT::AliGenCRT() : AliGenerator(-1)
82 // Set the origin above the vertex, on the surface.
102 fUnfoldedMomentumDist = 0;
107 //_____________________________________________________________________________
108 AliGenCRT::AliGenCRT(Int_t npart)
109 : AliGenerator(npart)
115 fTitle = "Cosmic Muons generator";
117 // Generate Muon- by default
126 // Set the origin above the vertex, on the surface.
128 fOrigin[1] = AliCRTConstants::fgDepth; // At the surface by default.
131 fZenithMin = 0.; // Generate veritcals by default.
137 fPResolution = 1.; // 1 GeV by default.
138 fPRange = 0.; // 0 GeV by default.
146 fUnfoldedMomentumDist = 0;
151 //_____________________________________________________________________________
152 AliGenCRT::AliGenCRT(const AliGenCRT& gen) : AliGenerator()
160 //_____________________________________________________________________________
161 AliGenCRT& AliGenCRT::operator= (const AliGenCRT& gen)
170 //_____________________________________________________________________________
171 AliGenCRT::~AliGenCRT()
176 if ( fAp ) delete fAp;
177 if ( fMomentumDist ) delete fMomentumDist;
178 if ( fUnfoldedMomentumDist ) delete fUnfoldedMomentumDist;
179 if ( fZenithDist ) delete fZenithDist;
180 if ( fPDist ) delete fPDist;
183 //_____________________________________________________________________________
184 void AliGenCRT::Generate()
187 // Generate on one trigger
190 Float_t polar[3]= {0,0,0}; // Polarization parameters
196 Float_t zenith, azimuth;
201 // Check if you need to use a distribution function for the momentum
202 if ( fMomentumDist ) {
203 pmom = - this->GetMomentum(); // Always downwards.
210 zenith = this->GetZenithAngle(pmom); // In degrees
212 // Take the azimuth random.
214 azimuth = (fAzimuthMin + (fAzimuthMax-fAzimuthMin)*random[0]); // In degrees
216 origin[0] = AliCRTConstants::fgDepth*TMath::Tan(zenith*kDegrad)*
217 TMath::Sin(azimuth*kDegrad);
218 origin[1] = AliCRTConstants::fgDepth;
219 origin[2] = AliCRTConstants::fgDepth*TMath::Tan(zenith*kDegrad)*
220 TMath::Cos(azimuth*kDegrad);
222 if ( fVertexSmear == kPerEvent ) {
226 // Don't smear the vertical position.
227 origin[j] = AliCRTConstants::fgDepth;
229 origin[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
230 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
236 if ( fCRMode == kSingleMuons ) {
239 // Generate kinematics a la AliGenBox
242 for(i=0;i<fNpart;i++) {
245 if(TestBit(kMomentumRange)) {
246 pmom = -( fPMin + random[1]*(fPMax - fPMin) ); // always downwards.
247 pt=pmom*TMath::Sin(zenith*kDegrad);
250 pt= -(fPtMin+random[1]*(fPtMax-fPtMin)); // always downwards.
251 pmom=pt/TMath::Sin(zenith*kDegrad);
254 p[0] = pt*TMath::Sin(azimuth*kDegrad);
255 p[1] = pmom*TMath::Cos(zenith*kDegrad);
256 p[2] = pt*TMath::Cos(azimuth*kDegrad);
258 if(fVertexSmear==kPerTrack) {
262 origin[j] = AliCRTConstants::fgDepth;
264 origin[j]=fOrigin[j]+fOsigma[j]*
265 TMath::Cos(2*random[2*j]*TMath::Pi())*
266 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
273 // End of generation a la AliGenBox
275 } else if ( fCRMode == kMuonBundle ) {
278 // Generate kinematics a la AliGenScan
294 for (Int_t ix=0; ix<fNx; ix++) {
295 for (Int_t iz=0; iz<fNz; iz++){
297 origin[0]+=ix*dx+2*(random[0]-0.5)*fOsigma[0];
298 origin[2]+=iz*dz+2*(random[1]-0.5)*fOsigma[2];
300 pmom = -(fPMin + random[3] *(fPMax - fPMax) ); // Always downward.
301 p[0] = TMath::Sin(zenith*kDegrad)*TMath::Sin(azimuth*kDegrad)*pmom;
302 p[1] = TMath::Cos(zenith*kDegrad)*pmom;
303 p[2] = TMath::Sin(zenith*kDegrad)*TMath::Cos(azimuth*kDegrad)*pmom;
308 // End of generation a la AliGenScan
310 } else if ( fCRMode == kMuonFlux ) {
312 // Always generate the things downwards.
313 p[0] = TMath::Sin(zenith*kDegrad)*TMath::Sin(azimuth*kDegrad)*pmom;
314 p[1] = TMath::Cos(zenith*kDegrad)*pmom;
315 p[2] = TMath::Sin(zenith*kDegrad)*TMath::Cos(azimuth*kDegrad)*pmom;
318 Error("Generate", "Generation Mode unknown!\n");
321 // Put the track on the stack.
322 SetTrack(fTrackIt,-1,fIpart,p,origin,polar,0,kPPrimary,nt);
327 //_____________________________________________________________________________
328 void AliGenCRT::Init()
331 // Initialize some internal methods.
334 // Determine some specific data members.
335 fPRange = TMath::Abs(fPMax-fPMin);
337 if ( fCRMode == kSingleMuons ) {
338 fCRModeName = new TString("Single Muons");
339 // Initialisation, check consistency of selected ranges
340 if(TestBit(kPtRange)&&TestBit(kMomentumRange))
341 Fatal("Init","You should not set the momentum range and the pt range!");
343 if((!TestBit(kPtRange))&&(!TestBit(kMomentumRange)))
344 Fatal("Init","You should set either the momentum or the pt range!");
346 } else if ( fCRMode == kMuonBundle ) {
347 fCRModeName = new TString("Muon Bundles");
349 } else if ( fCRMode == kMuonFlux ) {
350 fCRModeName = new TString("Muon Fluxes");
351 // Initialize the ditribution functions.
352 this->InitMomentumGeneration();
353 this->InitZenithalAngleGeneration();
356 Fatal("Generate", "Generation Mode unknown!\n");
362 //____________________________________________________________________________
363 void AliGenCRT::SetGridRange(Int_t nx,Float_t xwidth, Int_t nz, Float_t zwidth)
367 // This data shuold be used for Muon bundles generation.
374 // Print a message about the use, if the Mode has not been set, or
375 // it has to a different Mode.
376 if ( fCRMode != kMuonBundle ) {
377 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");
382 //____________________________________________________________________________
383 void AliGenCRT::InitApWeightFactors()
386 // This factors will be to correct the zenithal angle distribution
387 // acording the momentum
390 // Fill the array for the flux zenith angle dependence.
391 // at the index 0 of fAp[] will be the "factor" if we have a muon
393 Float_t a[6] = {-1.61, -1.50, -1.28, -0.94, -0.61, -0.22};
394 Float_t p[6] = { 0., 10., 30., 100., 300., 1000.};
396 // Get the information from the momentum
397 Int_t pEnd = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
398 // Initialize the Array of floats to hold the a(p) factors.
399 fAp = new TArrayF(pEnd);
403 for (Int_t i = 0; i < pEnd; i++ ) {
404 Float_t currentP = ((Float_t)i)*fPResolution;
405 if ( currentP < p[1] ) index = 0;
406 else if ( currentP >= p[1] && currentP < p[2] ) index = 1;
407 else if ( currentP >= p[2] && currentP < p[3] ) index = 2;
408 else if ( currentP >= p[3] && currentP < p[4] ) index = 3;
409 else if ( currentP >= p[4] ) index = 4;
411 Float_t ap = (currentP -p[index])*(a[index+1] - a[index])/
412 (p[index+1] - p[index]) + a[index];
418 //___________________________________________________________________________
419 void AliGenCRT::InitMomentumGeneration()
422 // Initialize a funtion to generate the momentum randomly
423 // acording this function.
426 // Check if we nned to initialize the function
427 if ( fPMin != fPMax ) {
429 // Check also if the function have been defined yet.
430 if ( !fMomentumDist ) {
432 // If not, use this function
433 const char* y = "log10(x)";
435 const char* h1Coef = "[0]*( %s*%s*%s/2 - (5*%s*%s/2) + 3*%s )";
436 const char* h2Coef = "[1]*( (-2*%s*%s*%s/3) + (3*%s*%s) - 10*%s/3 + 1 )";
437 const char* h3Coef = "[2]*( %s*%s*%s/6 - %s*%s/2 + %s/3 )";
438 const char* s2Coef = "[3]*( %s*%s*%s/3 - 2*%s*%s + 11*%s/3 - 2 )";
440 const char* h = "%s + %s + %s + %s";
441 const char* flux = "pow(10., %s)";
442 const char* normalizedFlux = "0.86*x*x*x*pow(10., %s)";
443 const char* paramNames[4] = {"H1", "H2", "H3", "S1"};
453 sprintf(buffer1, h1Coef, y, y, y, y, y, y);
454 sprintf(buffer2, h2Coef, y, y, y, y, y, y);
455 sprintf(buffer3, h3Coef, y, y, y, y, y, y);
456 sprintf(buffer4, s2Coef, y, y, y, y, y, y);
458 sprintf(buffer5, h, buffer1, buffer2, buffer3, buffer4);
460 sprintf(buffer6, flux, buffer5);
462 fMomentumDist = new TF1("fMomentumDist", buffer6, fPMin, fPMax);
463 sprintf(buffer7, normalizedFlux, buffer5);
464 fUnfoldedMomentumDist = new TF1("fUnfoldedMomentumDist", buffer7, fPMin, fPMax);
465 for (Int_t i = 0; i < 4; i++ ) {
466 fMomentumDist->SetParName(i, paramNames[i]);
467 fUnfoldedMomentumDist->SetParName(i, paramNames[i]);
470 fMomentumDist->SetParameter(0, 0.133);
471 fMomentumDist->SetParameter(1, -2.521);
472 fMomentumDist->SetParameter(2, -5.78);
473 fMomentumDist->SetParameter(3, -2.11);
475 fUnfoldedMomentumDist->SetParameter(0, 0.133);
476 fUnfoldedMomentumDist->SetParameter(1, -2.521);
477 fUnfoldedMomentumDist->SetParameter(2, -5.78);
478 fUnfoldedMomentumDist->SetParameter(3, -2.11);
486 //____________________________________________________________________________
487 void AliGenCRT::InitZenithalAngleGeneration()
490 // Initalize a distribution function for the zenith angle.
491 // This angle will be obtained randomly acording this function.
492 // The generated angles will been in degrees.
494 // Check if we need to create the function.
495 if ( fZenithMin != fZenithMax ) {
497 // Check also if another function have been defined.
498 if ( !fZenithDist ) {
500 // initialize the momentum dependent coefficients, a(p)
501 InitApWeightFactors();
503 // Define the standard function.
504 char* zenithalDisributionFunction = "1 + [0]*(1 - cos(x*3.14159265358979312/180))";
506 Int_t pEnd = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
507 fPDist = new TClonesArray("TF1", pEnd);
508 TClonesArray &angle = *fPDist;
509 for ( Int_t i = 0; i < pEnd; i++ ) {
510 // Fill the distribution
511 TF1* zenith = new(angle[i]) TF1("zenith",zenithalDisributionFunction, fZenithMin, fZenithMax);
513 // Fix the momentum dependent coefficients
514 zenith->SetParName(0, "a(p)");
515 zenith->SetParameter(0, fAp->At(i));
525 //____________________________________________________________________________
526 const Float_t AliGenCRT::GetZenithAngle(Float_t mom)
530 // Check if you need to generate a constant zenith angle.
531 if ( !fZenithDist ) {
532 // Check if you have defined an array of momentum functions
534 Int_t pIndex = TMath::Abs(TMath::Nint(mom));
535 TF1* zenithAngle = (TF1*)fPDist->UncheckedAt(pIndex);
536 zenith = zenithAngle->GetRandom();
540 if ( fCRMode != kMuonFlux ) {
541 // If you aren't generating muons obeying any ditribution
542 // only generate a flat zenith angle, acording the input settings
545 zenith = fZenithMin + (fZenithMax - fZenithMin)*random[0];
548 // Even if you are generating muons acording some distribution,
549 // but you don't want to ...
556 zenith = fZenithDist->GetRandom();
563 //____________________________________________________________________________