ACORDE/AliGenACORDE.cxx

   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 ClassImp(AliGenACORDE)
  50
  51 //_____________________________________________________________________________
  52 AliGenACORDE::AliGenACORDE()
  53   : AliGenerator(),
  54     fIpart(0),
  55     fCRMode(kSingleMuons),
  56     fCRModeName(0),
  57     fXwidth(0),
  58     fNx(1),
  59     fZwidth(0),
  60     fNz(1),
  61     fMuonGrid(kFALSE),
  62     fZenithMin(0),
  63     fZenithMax(0),
  64     fAzimuthMin(0),
  65     fAzimuthMax(0),
  66     fPRange(0),
  67     fPResolution(1),
  68     fAp(0),
  69     fMomentumDist(0),
  70     fUnfoldedMomentumDist(0),
  71     fZenithDist(0),
  72     fPDist(0),
  73     fNParticles(0)
  74 {
  75   //
  76   // Default ctor.
  77   //
  78 }
  79
  80 //_____________________________________________________________________________
  81 AliGenACORDE::AliGenACORDE(Int_t npart)
  82   : AliGenerator(npart),
  83     fIpart(kMuonMinus),
  84     fCRMode(kSingleMuons),
  85     fCRModeName(0),
  86     fXwidth(0),
  87     fNx(1),
  88     fZwidth(0),
  89     fNz(1),
  90     fMuonGrid(kFALSE),
  91     fZenithMin(0),
  92     fZenithMax(0),
  93     fAzimuthMin(0),
  94     fAzimuthMax(0),
  95     fPRange(0),
  96     fPResolution(1),
  97     fAp(0),
  98     fMomentumDist(0),
  99     fUnfoldedMomentumDist(0),
 100     fZenithDist(0),
 101     fPDist(0),
 102     fNParticles(0)
 103 {
 104   //
 105   // Standard ctor.
 106   //
 107   fName = "ACORDE";
 108   fTitle = "Cosmic Muons generator";
 109
 110   // Set the origin above the vertex, on the surface.
 111   fOrigin[0] = 0.;
 112   fOrigin[1] = AliACORDEConstants::Instance()->Depth(); // At the surface by default.
 113   fOrigin[2] = 0.;
 114 }
 115
 116 //_____________________________________________________________________________
 117 AliGenACORDE::~AliGenACORDE()
 118 {
 119   //
 120   // Default dtor.
 121   //
 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;
 127 }
 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   for (Int_t i = 0; i < fNParticles; i++ ) {
 138
 139     if ( fCRMode == kMuonBundle ) {
 140       this->GenerateOneMuonBundle();
 141
 142     } else if ( fCRMode == kSingleMuons ) {
 143       this->GenerateOneSingleMuon(kTRUE);
 144
 145     } else {
 146       // Generate only single muons following the parametrizations
 147       // for atmospheric muons.
 148       this->GenerateOneSingleMuon();
 149
 150     }
 151
 152   }
 153 }
 154
 155 //_____________________________________________________________________________
 156 void AliGenACORDE::Init()
 157 {
 158   //
 159   // Initialize some internal methods.
 160   //
 161
 162
 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");
 167
 168   // Determine some specific data members.
 169   fPRange = TMath::Abs(fPMax-fPMin);
 170
 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!");
 176
 177     if((!TestBit(kPtRange))&&(!TestBit(kMomentumRange)))
 178       Fatal("Init","You should set either the momentum or the pt range!");
 179
 180   } else if ( fCRMode == kMuonBundle ) {
 181     fCRModeName = new TString("Muon Bundles");
 182
 183   } else if ( fCRMode == kMuonFlux ) {
 184     fCRModeName = new TString("Muon Fluxes");
 185     // Initialize the ditribution functions.
 186     this->InitMomentumGeneration();
 187     this->InitZenithalAngleGeneration();
 188
 189   } else {
 190     Fatal("Generate", "Generation Mode unknown!\n");
 191
 192   }
 193
 194 }
 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   Float_t polar[3]= {0,0,0}; // Polarization parameters
 213   Float_t origin[3];
 214   Int_t nt;
 215   Float_t p[3];
 216   Float_t pmom, pt;
 217   Float_t random[6];
 218
 219   // Take the azimuth random.
 220   Rndm(random, 2);
 221   Float_t azimuth = fAzimuthMin + (fAzimuthMax-fAzimuthMin)*random[0];
 222   Float_t zenith = fZenithMin + (fZenithMax - fZenithMin)*random[1];
 223
 224   if ( withFlatMomentum ) {
 225     Rndm(random,3);
 226     if(TestBit(kMomentumRange)) {
 227       pmom = -( fPMin + random[0]*(fPMax - fPMin) ); // always downwards.
 228       pt = pmom*TMath::Sin(zenith*kDegrad);
 229     } else {
 230       pt = -( fPtMin + random[1]*(fPtMax - fPtMin)); // always downwards.
 231       pmom = pt/TMath::Sin(zenith*kDegrad);
 232     }
 233
 234   } else {
 235     if ( fMomentumDist ) {
 236       pmom = -this->GetMomentum(); // Always downwards.
 237     } else {
 238       pmom = -fPMin;
 239     }
 240     zenith = this->GetZenithAngle(pmom);  // In degrees
 241     pt = pmom*TMath::Sin(zenith*kDegrad);
 242   }
 243
 244   p[0] = pt*TMath::Sin(azimuth*kDegrad);
 245   p[1] = pmom*TMath::Cos(zenith*kDegrad);
 246   p[2] = pt*TMath::Cos(azimuth*kDegrad);
 247
 248   // Finaly the origin, with the smearing
 249   Rndm(random,6);
 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]));
 253
 254   origin[1] = AliACORDEConstants::Instance()->Depth();
 255
 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]));
 259
 260   // Put the track on the stack.
 261   PushTrack(fTrackIt,-1,fIpart,p,origin,polar,0,kPPrimary,nt);
 262
 263 }
 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   Float_t polar[3]= {0,0,0}; // Polarization parameters
 275   Float_t origin[3];
 276   Float_t p[3];
 277   Int_t nt;
 278   Float_t pmom;
 279   Float_t random[6];
 280
 281   Rndm(random, 3);
 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;
 286
 287   // Generate the kinematics a la AliGenScan (Andreas Morchs)
 288   Float_t dx, dz;
 289   if ( fNx > 0 ) {
 290     dx = fXwidth/fNx;
 291   } else {
 292     dx = 1e10;
 293     //dx = 100.;
 294   }
 295
 296   if ( fNz > 0 ) {
 297     dz = fZwidth/fNz;
 298   } else {
 299     dz = 1e10;
 300     //dz = 100.;
 301   }
 302
 303   origin[0] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
 304               TMath::Sin(azimuth*kDegrad);
 305   //origin[0] = 0.;
 306   origin[1] = AliACORDEConstants::Instance()->Depth();
 307   //origin[1] = 900;
 308   origin[2] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
 309               TMath::Cos(azimuth*kDegrad);
 310     //origin[2] = 0.;
 311
 312   for (Int_t ix = 0; ix < fNx; ix++ ) {
 313     for (Int_t iz = 0; iz < fNz; iz++ ) {
 314       Rndm(random,6);
 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];
 319       }
 320       if ( random[5] < 0.5 ) {
 321         origin[2] = -1*origin[2];
 322       }
 323
 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;
 328
 329       PushTrack(fTrackIt, -1, fIpart, p, origin, polar, 0, kPPrimary, nt);
 330     }
 331
 332   }
 333
 334 }
 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   fXwidth=xwidth;
 344   fNx=nx;
 345   fZwidth=zwidth;
 346   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   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");
 352     fMuonGrid = kTRUE;
 353   }
 354 }
 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   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.};
 369
 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);
 374
 375   Int_t index = 0;
 376
 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;
 384
 385     Float_t ap = (currentP -p[index])*(a[index+1] - a[index])/
 386                  (p[index+1] - p[index]) + a[index];
 387     fAp->AddAt(ap, i);
 388   }
 389
 390 }
 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   if ( fPMin != fPMax ) {
 402
 403     // Check also if the function have been defined yet.
 404     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       const char* paramNames[4] = {"H1", "H2", "H3", "S1"};
 418
 419       char buffer1[1024];
 420       char buffer2[1024];
 421       char buffer3[1024];
 422       char buffer4[1024];
 423       char buffer5[1024];
 424       char buffer6[1024];
 425       char buffer7[1024];
 426
 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);
 431
 432       snprintf(buffer5, 1023, h, buffer1, buffer2, buffer3, buffer4);
 433
 434       snprintf(buffer6, 1023, flux, buffer5);
 435
 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]);
 442       }
 443
 444       fMomentumDist->SetParameter(0, 0.133);
 445       fMomentumDist->SetParameter(1, -2.521);
 446       fMomentumDist->SetParameter(2, -5.78);
 447       fMomentumDist->SetParameter(3, -2.11);
 448
 449       fUnfoldedMomentumDist->SetParameter(0, 0.133);
 450       fUnfoldedMomentumDist->SetParameter(1, -2.521);
 451       fUnfoldedMomentumDist->SetParameter(2, -5.78);
 452       fUnfoldedMomentumDist->SetParameter(3, -2.11);
 453
 454     }
 455
 456   }
 457
 458 }
 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   if ( fZenithMin != fZenithMax ) {
 470
 471     // Check also if another function have been defined.
 472     if ( !fZenithDist ) {
 473
 474       // initialize the momentum dependent coefficients, a(p)
 475       this->InitApWeightFactors();
 476
 477       Int_t pEnd  = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
 478       char name[26];
 479       char title[52];
 480       fPDist = new TClonesArray("TH1F", pEnd);
 481       TClonesArray &mom = *fPDist;
 482       TH1F* zenith = 0;
 483       Float_t weight = 0;
 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));
 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         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);
 497         }
 498
 499       }
 500
 501     }
 502
 503   }
 504
 505 }
 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   if ( !fZenithDist ) {
 514     // Check if you have defined an array of momentum functions
 515     if ( fPDist ) {
 516       Int_t pIndex = TMath::Abs(TMath::Nint(mom));
 517       TH1F* cosZenithAngle = (TH1F*)fPDist->UncheckedAt(pIndex);
 518       Float_t tmpzenith = TMath::ACos(cosZenithAngle->GetRandom());
 519       // Correct the value
 520       zenith = kRaddeg*tmpzenith;
 521       return zenith;
 522     } else {
 523
 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
 527         Float_t random[2];
 528         Rndm(random, 2);
 529         zenith = fZenithMin + (fZenithMax - fZenithMin)*random[0];
 530
 531       } else {
 532         // Even if you are generating muons acording some distribution,
 533         // but you don't want to ...
 534         zenith = fZenithMin;
 535
 536       }
 537
 538     }
 539   } else {
 540     zenith = fZenithDist->GetRandom();
 541   }
 542
 543   return zenith;
 544 }
 545
 546 //_____________________________________________________________________________
 547 Float_t AliGenACORDE::GetMomentum() const
 548 {
 549   //
 550   //
 551   //
 552   return fMomentumDist->GetRandom();
 553 }