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