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