CRT/AliGenCRT.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/AliGenCRTClass.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 "AliGenCRT.h"
  39
  40 #include <TMCProcess.h>
  41 #include <TPDGCode.h>
  42 #include <TClonesArray.h>
  43 #include <TF1.h>
  44
  45 #include "AliRun.h"
  46 #include "AliConst.h"
  47
  48 ClassImp(AliGenCRT)
  49
  50 //_____________________________________________________________________________
  51 AliGenCRT::AliGenCRT()
  52   : AliGenerator(),
  53     fIpart(0),
  54     fCRMode(kSingleMuons),
  55     fCRModeName(0),
  56     fXwidth(0),
  57     fNx(1),
  58     fZwidth(0),
  59     fNz(1),
  60     fMuonGrid(kFALSE),
  61     fZenithMin(0),
  62     fZenithMax(0),
  63     fAzimuthMin(0),
  64     fAzimuthMax(0),
  65     fPRange(0),
  66     fPResolution(1),
  67     fAp(0),
  68     fMomentumDist(0),
  69     fUnfoldedMomentumDist(0),
  70     fZenithDist(0),
  71     fPDist(0)
  72 {
  73   //
  74   // Default ctor.
  75   //
  76 }
  77
  78 //_____________________________________________________________________________
  79 AliGenCRT::AliGenCRT(Int_t npart)
  80   : AliGenerator(npart),
  81     fIpart(kMuonMinus),
  82     fCRMode(kSingleMuons),
  83     fCRModeName(0),
  84     fXwidth(0),
  85     fNx(1),
  86     fZwidth(0),
  87     fNz(1),
  88     fMuonGrid(kFALSE),
  89     fZenithMin(0),
  90     fZenithMax(0),
  91     fAzimuthMin(0),
  92     fAzimuthMax(0),
  93     fPRange(0),
  94     fPResolution(1),
  95     fAp(0),
  96     fMomentumDist(0),
  97     fUnfoldedMomentumDist(0),
  98     fZenithDist(0),
  99     fPDist(0)
 100 {
 101   //
 102   // Standard ctor.
 103   //
 104   fName = "CRT";
 105   fTitle = "Cosmic Muons generator";
 106
 107   // Set the origin above the vertex, on the surface.
 108   fOrigin[0] = 0.;
 109   fOrigin[1] = AliCRTConstants::Instance()->Depth(); // At the surface by default.
 110   fOrigin[2] = 0.;
 111 }
 112
 113 //_____________________________________________________________________________
 114 AliGenCRT::AliGenCRT(const AliGenCRT& gen)
 115   : AliGenerator(gen)
 116 {
 117   //
 118   // Copy constructor
 119   //
 120   gen.Copy(*this);
 121 }
 122
 123 //_____________________________________________________________________________
 124 AliGenCRT& AliGenCRT::operator=(const AliGenCRT& gen)
 125 {
 126   //
 127   // Asingment operator
 128   //
 129   gen.Copy(*this);
 130   return *this;
 131 }
 132
 133 //_____________________________________________________________________________
 134 AliGenCRT::~AliGenCRT()
 135 {
 136   //
 137   // Default dtor.
 138   //
 139   if ( fPDist ) {fPDist->Delete(); delete fPDist; fPDist = 0;}
 140   if ( fUnfoldedMomentumDist ) delete fUnfoldedMomentumDist;
 141   if ( fMomentumDist ) delete fMomentumDist;
 142   if ( fAp )           delete fAp;
 143   if ( fCRModeName )   delete fCRModeName;
 144 }
 145
 146 //_____________________________________________________________________________
 147 void AliGenCRT::Generate()
 148 {
 149   //
 150   // Generate on one trigger
 151   // Call the respective method inside the loop for the number
 152   // of tracks per trigger.
 153
 154   for (Int_t i = 0; i < fNpart; i++ ) {
 155
 156     if ( fCRMode == kMuonBundle ) {
 157       this->GenerateOneMuonBundle();
 158
 159     } else if ( fCRMode == kSingleMuons ) {
 160       this->GenerateOneSingleMuon(kTRUE);
 161
 162     } else {
 163       // Generate only single muons but following the parametrizations
 164       // for atmospheric muons.
 165       this->GenerateOneSingleMuon();
 166
 167     }
 168
 169   }
 170 }
 171
 172 //_____________________________________________________________________________
 173 void AliGenCRT::Init()
 174 {
 175   //
 176   // Initialize some internal methods.
 177   //
 178
 179   // Determine some specific data members.
 180   fPRange = TMath::Abs(fPMax-fPMin);
 181
 182   if ( fCRMode == kSingleMuons ) {
 183     fCRModeName = new TString("Single Muons");
 184     // Initialisation, check consistency of selected ranges
 185     if(TestBit(kPtRange)&&TestBit(kMomentumRange))
 186       Fatal("Init","You should not set the momentum range and the pt range!");
 187
 188     if((!TestBit(kPtRange))&&(!TestBit(kMomentumRange)))
 189       Fatal("Init","You should set either the momentum or the pt range!");
 190
 191   } else if ( fCRMode == kMuonBundle ) {
 192     fCRModeName = new TString("Muon Bundles");
 193
 194   } else if ( fCRMode == kMuonFlux ) {
 195     fCRModeName = new TString("Muon Fluxes");
 196     // Initialize the ditribution functions.
 197     this->InitMomentumGeneration();
 198     this->InitZenithalAngleGeneration();
 199
 200   } else {
 201     Fatal("Generate", "Generation Mode unknown!\n");
 202
 203   }
 204
 205 }
 206
 207 //____________________________________________________________________________
 208 void AliGenCRT::GenerateOneSingleMuon(Bool_t withFlatMomentum)
 209 {
 210   //
 211   // Generate One Single Muon
 212   // This method will generate only one muon.
 213   // The momentum will be randomly flat distributed if
 214   // the paremeter "withFlatMomentum" is set to kTRUE,
 215   // otherwise the momentum will generate acordingly the parametrization
 216   // given by
 217   // and adpted from Bruno Alessandro's implementation with the
 218   // CERNLIB to AliRoot.
 219   // The "withFlatMomentum" parameter also will be used to generate
 220   // the muons with a flat Zenithal angle distribution.
 221   // Do the smearing here, so that means per track.
 222
 223   Float_t polar[3]= {0,0,0}; // Polarization parameters
 224   Float_t origin[3];
 225   Int_t nt;
 226   Float_t p[3];
 227   Float_t pmom, pt;
 228   Float_t random[6];
 229
 230   // Take the azimuth random.
 231   Rndm(random, 2);
 232   Float_t azimuth = fAzimuthMin + (fAzimuthMax-fAzimuthMin)*random[0];
 233   Float_t zenith = fZenithMin + (fZenithMax - fZenithMin)*random[1];
 234
 235   if ( withFlatMomentum ) {
 236     Rndm(random,3);
 237     if(TestBit(kMomentumRange)) {
 238       pmom = -( fPMin + random[0]*(fPMax - fPMin) ); // always downwards.
 239       pt = pmom*TMath::Sin(zenith*kDegrad);
 240     } else {
 241       pt = -( fPtMin + random[1]*(fPtMax - fPtMin)); // always downwards.
 242       pmom = pt/TMath::Sin(zenith*kDegrad);
 243     }
 244
 245   } else {
 246     if ( fMomentumDist ) {
 247       pmom = -this->GetMomentum(); // Always downwards.
 248     } else {
 249       pmom = -fPMin;
 250     }
 251     zenith = this->GetZenithAngle(pmom);  // In degrees
 252     pt = pmom*TMath::Sin(zenith*kDegrad);
 253   }
 254
 255   p[0] = pt*TMath::Sin(azimuth*kDegrad);
 256   p[1] = pmom*TMath::Cos(zenith*kDegrad);
 257   p[2] = pt*TMath::Cos(azimuth*kDegrad);
 258
 259   // Finaly the origin, with the smearing
 260   Rndm(random,6);
 261   origin[0] = AliCRTConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
 262     TMath::Sin(azimuth*kDegrad);
 263   // + fOsigma[0]* TMath::Cos(2*random[0]*TMath::Pi())*TMath::Sqrt(-2*TMath::Log(random[1]));
 264
 265   origin[1] = AliCRTConstants::Instance()->Depth();
 266
 267   origin[2] = AliCRTConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
 268     TMath::Cos(azimuth*kDegrad);
 269   // + fOsigma[2]* TMath::Cos(2*random[2]*TMath::Pi())*TMath::Sqrt(-2*TMath::Log(random[3]));;
 270
 271   // Put the track on the stack.
 272   PushTrack(fTrackIt,-1,fIpart,p,origin,polar,0,kPPrimary,nt);
 273
 274 }
 275
 276 //____________________________________________________________________________
 277 void AliGenCRT::GenerateOneMuonBundle()
 278 {
 279   //
 280   // Generate One Muon Bundle method
 281   // This method will generate a bunch of muons following the
 282   // procedure of the AliGenScan class.
 283   // These muons will be generated with flat momentum.
 284
 285   Float_t polar[3]= {0,0,0}; // Polarization parameters
 286   Float_t origin[3];
 287   Float_t p[3];
 288   Int_t nt;
 289   Float_t pmom;
 290   Float_t random[6];
 291
 292   Rndm(random, 3);
 293   //Float_t zenith = fZenithMin + (fZenithMax - fZenithMin)*random[1];
 294   //Float_t azimuth = fAzimuthMin + (fAzimuthMax-fAzimuthMin)*random[2];
 295   Float_t zenith = 10;
 296   Float_t azimuth = 30;
 297
 298   // Generate the kinematics a la AliGenScan (Andreas Morchs)
 299   Float_t dx, dz;
 300   if ( fNx > 0 ) {
 301     dx = fXwidth/fNx;
 302   } else {
 303     dx = 1e10;
 304     //dx = 100.;
 305   }
 306
 307   if ( fNz > 0 ) {
 308     dz = fZwidth/fNz;
 309   } else {
 310     dz = 1e10;
 311     //dz = 100.;
 312   }
 313
 314   //origin[0] = AliCRTConstants::fgDepth*TMath::Tan(zenith*kDegrad)*
 315   origin[0] = 0.;
 316     //TMath::Sin(azimuth*kDegrad);
 317   //origin[1] = AliCRTConstants::fgDepth;
 318   origin[1] = 900;
 319   //origin[2] = AliCRTConstants::fgDepth*TMath::Tan(zenith*kDegrad)*
 320   origin[2] = 0.;
 321     //TMath::Cos(azimuth*kDegrad);
 322
 323   for (Int_t ix = 0; ix < fNx; ix++ ) {
 324     for (Int_t iz = 0; iz < fNz; iz++ ) {
 325       Rndm(random,6);
 326       origin[0]+=ix*dx+2*(random[1]-0.5)*fOsigma[0];
 327       origin[2]+=iz*dz+2*(random[2]-0.5)*fOsigma[2];
 328       if ( random[4] < 0.5 ) {
 329         origin[0] = -1*origin[0];
 330       }
 331       if ( random[5] < 0.5 ) {
 332         origin[2] = -1*origin[2];
 333       }
 334
 335       pmom = -(fPMin + random[3] *(fPMax - fPMax) ); // Always downwards
 336       p[0] = TMath::Sin(zenith*kDegrad)*TMath::Sin(azimuth*kDegrad)*pmom;
 337       p[1] = TMath::Cos(zenith*kDegrad)*pmom;
 338       p[2] = TMath::Sin(zenith*kDegrad)*TMath::Cos(azimuth*kDegrad)*pmom;
 339
 340       PushTrack(fTrackIt, -1, fIpart, p, origin, polar, 0, kPPrimary, nt);
 341     }
 342
 343   }
 344
 345 }
 346
 347 //____________________________________________________________________________
 348 void AliGenCRT::SetGridRange(Int_t nx,Float_t xwidth, Int_t nz, Float_t zwidth)
 349 {
 350   //
 351   // Define the grid
 352   // This data shuold be used for Muon bundles generation.
 353   //
 354   fXwidth=xwidth;
 355   fNx=nx;
 356   fZwidth=zwidth;
 357   fNz=nz;
 358
 359   // Print a message  about the use, if the Mode has not been set, or
 360   // it has to a different Mode.
 361   if ( fCRMode != kMuonBundle ) {
 362     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");
 363     fMuonGrid = kTRUE;
 364   }
 365 }
 366
 367 //____________________________________________________________________________
 368 void AliGenCRT::InitApWeightFactors()
 369 {
 370   //
 371   // This factors will be  to correct the zenithal angle distribution
 372   // acording the momentum
 373
 374   //
 375   // Fill the array for the flux zenith angle dependence.
 376   // at the index 0 of fAp[] will be the "factor" if we have a muon
 377   // of 0 GeV.
 378   Float_t a[6] = {-1.61, -1.50, -1.28, -0.94, -0.61, -0.22};
 379   Float_t p[6] = { 0., 10., 30., 100., 300., 1000.};
 380
 381   // Get the information from the momentum
 382   Int_t pEnd  = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
 383   // Initialize the Array of floats to hold the a(p) factors.
 384   fAp = new TArrayF(pEnd);
 385
 386   Int_t index = 0;
 387
 388   for (Int_t i = 0; i < pEnd; i++ ) {
 389     Float_t currentP = ((Float_t)i)*fPResolution;
 390     if ( currentP < p[1] )                          index = 0;
 391     else if ( currentP >= p[1] && currentP < p[2] ) index = 1;
 392     else if ( currentP >= p[2] && currentP < p[3] ) index = 2;
 393     else if ( currentP >= p[3] && currentP < p[4] ) index = 3;
 394     else if ( currentP >= p[4] )                    index = 4;
 395
 396     Float_t ap = (currentP -p[index])*(a[index+1] - a[index])/
 397                  (p[index+1] - p[index]) + a[index];
 398     fAp->AddAt(ap, i);
 399   }
 400
 401 }
 402
 403 //___________________________________________________________________________
 404 void AliGenCRT::InitMomentumGeneration()
 405 {
 406   //
 407   // Initialize a funtion to generate the momentum randomly
 408   // acording this function.
 409   //
 410
 411   // Check if we nned to initialize the function
 412   if ( fPMin != fPMax ) {
 413
 414     // Check also if the function have been defined yet.
 415     if ( !fMomentumDist ) {
 416
 417       // If not, use this function
 418       const char* y      = "log10(x)";
 419
 420       const char* h1Coef = "[0]*( %s*%s*%s/2 - (5*%s*%s/2) + 3*%s )";
 421       const char* h2Coef = "[1]*( (-2*%s*%s*%s/3) + (3*%s*%s) - 10*%s/3 + 1 )";
 422       const char* h3Coef = "[2]*( %s*%s*%s/6 - %s*%s/2 + %s/3 )";
 423       const char* s2Coef = "[3]*( %s*%s*%s/3 - 2*%s*%s + 11*%s/3 - 2 )";
 424
 425       const char* h = "%s + %s + %s + %s";
 426       const char* flux = "pow(10., %s)";
 427       const char* normalizedFlux = "0.86*x*x*x*pow(10., %s)";
 428       const char* paramNames[4] = {"H1", "H2", "H3", "S1"};
 429
 430       char buffer1[1024];
 431       char buffer2[1024];
 432       char buffer3[1024];
 433       char buffer4[1024];
 434       char buffer5[1024];
 435       char buffer6[1024];
 436       char buffer7[1024];
 437
 438       sprintf(buffer1, h1Coef, y, y, y, y, y, y);
 439       sprintf(buffer2, h2Coef, y, y, y, y, y, y);
 440       sprintf(buffer3, h3Coef, y, y, y, y, y, y);
 441       sprintf(buffer4, s2Coef, y, y, y, y, y, y);
 442
 443       sprintf(buffer5, h, buffer1, buffer2, buffer3, buffer4);
 444
 445       sprintf(buffer6, flux, buffer5);
 446
 447       fMomentumDist = new TF1("fMomentumDist", buffer6, fPMin, fPMax);
 448       sprintf(buffer7, normalizedFlux, buffer5);
 449       fUnfoldedMomentumDist = new TF1("fUnfoldedMomentumDist", buffer7, fPMin, fPMax);
 450       for (Int_t i = 0; i < 4; i++ ) {
 451         fMomentumDist->SetParName(i, paramNames[i]);
 452         fUnfoldedMomentumDist->SetParName(i, paramNames[i]);
 453       }
 454
 455       fMomentumDist->SetParameter(0, 0.133);
 456       fMomentumDist->SetParameter(1, -2.521);
 457       fMomentumDist->SetParameter(2, -5.78);
 458       fMomentumDist->SetParameter(3, -2.11);
 459
 460       fUnfoldedMomentumDist->SetParameter(0, 0.133);
 461       fUnfoldedMomentumDist->SetParameter(1, -2.521);
 462       fUnfoldedMomentumDist->SetParameter(2, -5.78);
 463       fUnfoldedMomentumDist->SetParameter(3, -2.11);
 464
 465     }
 466
 467   }
 468
 469 }
 470
 471 //____________________________________________________________________________
 472 void AliGenCRT::InitZenithalAngleGeneration()
 473 {
 474   //
 475   // Initalize a distribution function for the zenith angle.
 476   // This angle will be obtained randomly acording this function.
 477   // The generated angles  will been in degrees.
 478
 479   // Check if we need to create the function.
 480   if ( fZenithMin != fZenithMax ) {
 481
 482     // Check also if another function have been defined.
 483     if ( !fZenithDist ) {
 484
 485       // initialize the momentum dependent coefficients, a(p)
 486       this->InitApWeightFactors();
 487
 488       // Define the standard function.
 489       const char* zenithalDisributionFunction = "1 + [0]*(1 - cos(x*3.14159265358979312/180))";
 490
 491       Int_t pEnd  = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
 492       fPDist = new TClonesArray("TF1", pEnd);
 493       TClonesArray &angle = *fPDist;
 494       for ( Int_t i = 0; i < pEnd; i++ ) {
 495         // Fill the distribution
 496         TF1* zenith = new(angle[i]) TF1("zenith",zenithalDisributionFunction, fZenithMin, fZenithMax);
 497
 498         // Fix the momentum dependent coefficients
 499         zenith->SetParName(0, "a(p)");
 500         zenith->SetParameter(0, fAp->At(i));
 501
 502       }
 503
 504     }
 505
 506   }
 507
 508 }
 509
 510 //____________________________________________________________________________
 511 const Float_t AliGenCRT::GetZenithAngle(Float_t mom) const
 512 {
 513
 514   Float_t zenith = 0.;
 515   // Check if you need to generate a constant zenith angle.
 516   if ( !fZenithDist ) {
 517     // Check if you have defined an array of momentum functions
 518     if ( fPDist ) {
 519       Int_t pIndex = TMath::Abs(TMath::Nint(mom));
 520       TF1* zenithAngle = (TF1*)fPDist->UncheckedAt(pIndex);
 521
 522       Float_t tmpzenith = TMath::Cos(zenithAngle->GetRandom()*kDegrad);
 523       // Correct the value
 524       zenith = kRaddeg*TMath::ACos(1 - (tmpzenith - 1)/fAp->At(pIndex));
 525       return zenith;
 526     } else {
 527
 528       if ( fCRMode != kMuonFlux ) {
 529         // If you aren't generating muons obeying any ditribution
 530         // only generate a flat zenith angle, acording the input settings
 531         Float_t random[2];
 532         Rndm(random, 2);
 533         zenith = fZenithMin + (fZenithMax - fZenithMin)*random[0];
 534
 535       } else {
 536         // Even if you are generating muons acording some distribution,
 537         // but you don't want to ...
 538         zenith = fZenithMin;
 539
 540       }
 541
 542     }
 543   } else {
 544     zenith = fZenithDist->GetRandom();
 545   }
 546
 547   return zenith;
 548 }
 549
 550 //_____________________________________________________________________________
 551 const Float_t AliGenCRT::GetMomentum() const
 552 {
 553   //
 554   //
 555   //
 556   return fMomentumDist->GetRandom();
 557 }