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