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 /*
  17 $Log$
  18 Revision 1.4  2002/11/21 19:34:52  alibrary
  19 Removing AliMC and AliMCProcess
  20
  21 Revision 1.3  2002/10/23 06:47:56  alibrary
  22 Introducing Riostream.h
  23
  24 Revision 1.2  2002/10/14 14:55:34  hristov
  25 Merging the VirtualMC branch to the main development branch (HEAD)
  26
  27 Revision 1.1.2.1  2002/10/10 14:40:31  hristov
  28 Updating VirtualMC to v3-09-02
  29
  30 Revision 1.1  2002/10/07 11:25:28  gamez
  31 First version, generation of cosmic muons on the surface
  32
  33
  34 */
  35
  36 /////////////////////////////////////////////////////////////////////////////
  37 //
  38 //  Contain parametrizations to generate atmospheric muons, and also
  39 //  to generate single muons and muon bundles at surface level.
  40 //
  41 //Begin_Html
  42 /*
  43 <img src="picts/AliGenCRTClass.gif">
  44 </pre>
  45 <br clear=left>
  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>.
  49 </font>
  50 <pre>
  51 */
  52 //End_Html
  53 //
  54 /////////////////////////////////////////////////////////////////////////////
  55
  56 #include <Riostream.h>
  57 #include <TMCProcess.h>
  58 #include <TPDGCode.h>
  59
  60
  61 #include "AliConst.h"
  62 #include "AliGenCRT.h"
  63 #include "AliRun.h"
  64
  65 ClassImp(AliGenCRT)
  66
  67 //_____________________________________________________________________________
  68 AliGenCRT::AliGenCRT() : AliGenerator(-1)
  69 {
  70   //
  71   // Default ctor.
  72   //
  73   fIpart = 0;
  74
  75   fXwidth=0;
  76   fNx=1;
  77   fZwidth=0;
  78   fNz=1;
  79   fMuonGrid = kFALSE;
  80
  81
  82   // Set the origin above the vertex, on the surface.
  83   fOrigin[0] = 0.;
  84   fOrigin[1] = 0.;
  85   fOrigin[2] = 0.;
  86
  87   fZenithMin = 0.;
  88   fZenithMax = 0.;
  89
  90   fAzimuthMin = 0.;
  91   fAzimuthMax = 0.;
  92
  93   fPResolution = 1.;
  94   fPRange = 0.;
  95
  96   fMomentumDist = 0;
  97   fZenithDist = 0;
  98
  99   fPDist = 0;
 100
 101   fAp = 0;
 102   fUnfoldedMomentumDist = 0;
 103
 104   fCRModeName = 0;
 105 }
 106
 107 //_____________________________________________________________________________
 108 AliGenCRT::AliGenCRT(Int_t npart)
 109   : AliGenerator(npart)
 110 {
 111   //
 112   // Standard ctor.
 113   //
 114   fName = "CRT";
 115   fTitle = "Cosmic Muons generator";
 116
 117   // Generate Muon- by default
 118   fIpart = kMuonMinus;
 119
 120   fXwidth=0;
 121   fNx=1;
 122   fZwidth=0;
 123   fNz=1;
 124   fMuonGrid = kFALSE;
 125
 126   // Set the origin above the vertex, on the surface.
 127   fOrigin[0] = 0.;
 128   fOrigin[1] = AliCRTConstants::fgDepth; // At the surface by default.
 129   fOrigin[2] = 0.;
 130
 131   fZenithMin = 0.; // Generate veritcals by default.
 132   fZenithMax = 0.;
 133
 134   fAzimuthMin = 0.;
 135   fAzimuthMax = 0.;
 136
 137   fPResolution = 1.; // 1 GeV by default.
 138   fPRange = 0.;      // 0 GeV by default.
 139
 140   fMomentumDist = 0;
 141   fZenithDist = 0;
 142
 143   fPDist = 0;
 144
 145   fAp = 0;
 146   fUnfoldedMomentumDist = 0;
 147
 148   fCRModeName = 0;
 149 }
 150
 151 //_____________________________________________________________________________
 152 AliGenCRT::AliGenCRT(const AliGenCRT& gen) : AliGenerator()
 153 {
 154   //
 155   // Copy ctor.
 156   //
 157   gen.Copy(*this);
 158 }
 159
 160 //_____________________________________________________________________________
 161 AliGenCRT& AliGenCRT::operator= (const AliGenCRT& gen)
 162 {
 163   //
 164   // Asingment ctor.
 165   //
 166   gen.Copy(*this);
 167   return *this;
 168 }
 169
 170 //_____________________________________________________________________________
 171 AliGenCRT::~AliGenCRT()
 172 {
 173   //
 174   // Default dtor.
 175   //
 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;
 181 }
 182
 183 //_____________________________________________________________________________
 184 void AliGenCRT::Generate()
 185 {
 186   //
 187   // Generate on one trigger
 188   //
 189
 190   Float_t polar[3]= {0,0,0}; // Polarization parameters
 191   //
 192   Float_t origin[3];
 193   Float_t p[3];
 194   Int_t i, j, nt;
 195   Float_t pmom, pt;
 196   Float_t zenith, azimuth;
 197   //
 198   Float_t random[6];
 199
 200
 201   // Check if you need to use a distribution function for the momentum
 202   if ( fMomentumDist ) {
 203     pmom = - this->GetMomentum(); // Always downwards.
 204
 205   } else {
 206     pmom = -fPMin;
 207
 208   }
 209
 210   zenith = this->GetZenithAngle(pmom);  // In degrees
 211
 212   // Take the azimuth random.
 213   Rndm(random, 2);
 214   azimuth = (fAzimuthMin + (fAzimuthMax-fAzimuthMin)*random[0]); // In degrees
 215
 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);
 221
 222   if ( fVertexSmear == kPerEvent ) {
 223     Rndm(random,6);
 224     for (j=0;j<3;j++) {
 225       if ( j == 1 ) {
 226         // Don't smear the vertical position.
 227         origin[j] = AliCRTConstants::fgDepth;
 228       } else {
 229         origin[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
 230           TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
 231       }
 232     }
 233   }
 234
 235
 236   if ( fCRMode == kSingleMuons ) {
 237
 238     //
 239     // Generate kinematics a la AliGenBox
 240     //
 241
 242     for(i=0;i<fNpart;i++) {
 243       Rndm(random,3);
 244
 245       if(TestBit(kMomentumRange)) {
 246         pmom = -( fPMin + random[1]*(fPMax - fPMin) ); // always downwards.
 247         pt=pmom*TMath::Sin(zenith*kDegrad);
 248       } else {
 249
 250         pt= -(fPtMin+random[1]*(fPtMax-fPtMin)); // always downwards.
 251         pmom=pt/TMath::Sin(zenith*kDegrad);
 252       }
 253
 254       p[0] = pt*TMath::Sin(azimuth*kDegrad);
 255       p[1] = pmom*TMath::Cos(zenith*kDegrad);
 256       p[2] = pt*TMath::Cos(azimuth*kDegrad);
 257
 258       if(fVertexSmear==kPerTrack) {
 259         Rndm(random,6);
 260         for (j=0;j<3;j++) {
 261           if ( j == 1 ) {
 262             origin[j] = AliCRTConstants::fgDepth;
 263           } else {
 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]));
 267           }
 268         }
 269       }
 270
 271     }
 272
 273     // End of generation a la AliGenBox
 274
 275   } else if ( fCRMode == kMuonBundle ) {
 276
 277     //
 278     // Generate kinematics a la AliGenScan
 279     //
 280     Float_t dx,dz;
 281     //
 282     if (fNx > 0) {
 283       dx=fXwidth/fNx;
 284     } else {
 285       dx=1e10;
 286     }
 287
 288     if (fNz > 0) {
 289       dz=fZwidth/fNz;
 290     } else {
 291       dz=1e10;
 292     }
 293
 294     for (Int_t ix=0; ix<fNx; ix++) {
 295       for (Int_t iz=0; iz<fNz; iz++){
 296         Rndm(random,6);
 297         origin[0]+=ix*dx+2*(random[0]-0.5)*fOsigma[0];
 298         origin[2]+=iz*dz+2*(random[1]-0.5)*fOsigma[2];
 299
 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;
 304
 305       }
 306     }
 307
 308     // End of generation a la AliGenScan
 309
 310   } else if ( fCRMode == kMuonFlux ) {
 311
 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;
 316
 317   } else {
 318     Error("Generate", "Generation Mode unknown!\n");
 319   }
 320
 321   // Put the track on the stack.
 322   SetTrack(fTrackIt,-1,fIpart,p,origin,polar,0,kPPrimary,nt);
 323
 324
 325 }
 326
 327 //_____________________________________________________________________________
 328 void AliGenCRT::Init()
 329 {
 330   //
 331   // Initialize some internal methods.
 332   //
 333
 334   // Determine some specific data members.
 335   fPRange = TMath::Abs(fPMax-fPMin);
 336
 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!");
 342
 343     if((!TestBit(kPtRange))&&(!TestBit(kMomentumRange)))
 344       Fatal("Init","You should set either the momentum or the pt range!");
 345
 346   } else if ( fCRMode == kMuonBundle ) {
 347     fCRModeName = new TString("Muon Bundles");
 348
 349   } else if ( fCRMode == kMuonFlux ) {
 350     fCRModeName = new TString("Muon Fluxes");
 351     // Initialize the ditribution functions.
 352     this->InitMomentumGeneration();
 353     this->InitZenithalAngleGeneration();
 354
 355   } else {
 356     Fatal("Generate", "Generation Mode unknown!\n");
 357
 358   }
 359
 360 }
 361
 362 //____________________________________________________________________________
 363 void AliGenCRT::SetGridRange(Int_t nx,Float_t xwidth, Int_t nz, Float_t zwidth)
 364 {
 365   //
 366   // Define the grid
 367   // This data shuold be used for Muon bundles generation.
 368   //
 369   fXwidth=xwidth;
 370   fNx=nx;
 371   fZwidth=zwidth;
 372   fNz=nz;
 373
 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");
 378     fMuonGrid = kTRUE;
 379   }
 380 }
 381
 382 //____________________________________________________________________________
 383 void AliGenCRT::InitApWeightFactors()
 384 {
 385   //
 386   // This factors will be  to correct the zenithal angle distribution
 387   // acording the momentum
 388
 389   //
 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
 392   // of 0 GeV.
 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.};
 395
 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);
 400
 401   Int_t index = 0;
 402
 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;
 410
 411     Float_t ap = (currentP -p[index])*(a[index+1] - a[index])/
 412                  (p[index+1] - p[index]) + a[index];
 413     fAp->AddAt(ap, i);
 414   }
 415
 416 }
 417
 418 //___________________________________________________________________________
 419 void AliGenCRT::InitMomentumGeneration()
 420 {
 421   //
 422   // Initialize a funtion to generate the momentum randomly
 423   // acording this function.
 424   //
 425
 426   // Check if we nned to initialize the function
 427   if ( fPMin != fPMax ) {
 428
 429     // Check also if the function have been defined yet.
 430     if ( !fMomentumDist ) {
 431
 432       // If not, use this function
 433       const char* y      = "log10(x)";
 434
 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 )";
 439
 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"};
 444
 445       char buffer1[1024];
 446       char buffer2[1024];
 447       char buffer3[1024];
 448       char buffer4[1024];
 449       char buffer5[1024];
 450       char buffer6[1024];
 451       char buffer7[1024];
 452
 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);
 457
 458       sprintf(buffer5, h, buffer1, buffer2, buffer3, buffer4);
 459
 460       sprintf(buffer6, flux, buffer5);
 461
 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]);
 468       }
 469
 470       fMomentumDist->SetParameter(0, 0.133);
 471       fMomentumDist->SetParameter(1, -2.521);
 472       fMomentumDist->SetParameter(2, -5.78);
 473       fMomentumDist->SetParameter(3, -2.11);
 474
 475       fUnfoldedMomentumDist->SetParameter(0, 0.133);
 476       fUnfoldedMomentumDist->SetParameter(1, -2.521);
 477       fUnfoldedMomentumDist->SetParameter(2, -5.78);
 478       fUnfoldedMomentumDist->SetParameter(3, -2.11);
 479
 480     }
 481
 482   }
 483
 484 }
 485
 486 //____________________________________________________________________________
 487 void AliGenCRT::InitZenithalAngleGeneration()
 488 {
 489   //
 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.
 493
 494   // Check if we need to create the function.
 495   if ( fZenithMin != fZenithMax ) {
 496
 497     // Check also if another function have been defined.
 498     if ( !fZenithDist ) {
 499
 500       // initialize the momentum dependent coefficients, a(p)
 501       InitApWeightFactors();
 502
 503       // Define the standard function.
 504       char* zenithalDisributionFunction = "1 + [0]*(1 - cos(x*3.14159265358979312/180))";
 505
 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);
 512
 513         // Fix the momentum dependent coefficients
 514         zenith->SetParName(0, "a(p)");
 515         zenith->SetParameter(0, fAp->At(i));
 516
 517       }
 518
 519     }
 520
 521   }
 522
 523 }
 524
 525 //____________________________________________________________________________
 526 const Float_t AliGenCRT::GetZenithAngle(Float_t mom)
 527 {
 528
 529   Float_t zenith = 0.;
 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
 533     if ( fPDist ) {
 534       Int_t pIndex = TMath::Abs(TMath::Nint(mom));
 535       TF1* zenithAngle = (TF1*)fPDist->UncheckedAt(pIndex);
 536       zenith = zenithAngle->GetRandom();
 537       return zenith;
 538     } else {
 539
 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
 543         Float_t random[2];
 544         Rndm(random, 2);
 545         zenith = fZenithMin + (fZenithMax - fZenithMin)*random[0];
 546
 547       } else {
 548         // Even if you are generating muons acording some distribution,
 549         // but you don't want to ...
 550         zenith = fZenithMin;
 551
 552       }
 553
 554     }
 555   } else {
 556     zenith = fZenithDist->GetRandom();
 557   }
 558
 559   return zenith;
 560
 561 }
 562
 563 //____________________________________________________________________________