[u/mrichter/AliRoot.git] / ACORDE / AliGenACORDE.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

/////////////////////////////////////////////////////////////////////////////
//
//  Contain parametrizations to generate atmospheric muons, and also
//  to generate single muons and muon bundles at surface level.
//
//Begin_Html
/*
<img src="picts/AliGenACORDEClass.gif">
</pre>
<br clear=left>
<font size=+2 color=red>
<p>The responsible person for this module is
<a href="mailto:Enrique.Gamez.Flores@cern.ch">Enrique Gamez</a>.
</font>
<pre>
*/
//End_Html
//
/////////////////////////////////////////////////////////////////////////////

#include "AliGenACORDE.h"

#include <TMCProcess.h>
#include <TPDGCode.h>
#include <TClonesArray.h>
#include <TF1.h>
#include <TH1F.h>

#include "AliRun.h"
#include "AliConst.h"

ClassImp(AliGenACORDE)

//_____________________________________________________________________________
AliGenACORDE::AliGenACORDE()
  : AliGenerator(),
    fIpart(0),
    fCRMode(kSingleMuons),
    fCRModeName(0),
    fXwidth(0),
    fNx(1),
    fZwidth(0),
    fNz(1),
    fMuonGrid(kFALSE),
    fZenithMin(0),
    fZenithMax(0),
    fAzimuthMin(0),
    fAzimuthMax(0),
    fPRange(0),
    fPResolution(1),
    fAp(0),
    fMomentumDist(0),
    fUnfoldedMomentumDist(0),
    fZenithDist(0),
    fPDist(0)
{
  //
  // Default ctor.
  //
}

//_____________________________________________________________________________
AliGenACORDE::AliGenACORDE(Int_t npart) 
  : AliGenerator(npart),
    fIpart(kMuonMinus),
    fCRMode(kSingleMuons),
    fCRModeName(0),
    fXwidth(0),
    fNx(1),
    fZwidth(0),
    fNz(1),
    fMuonGrid(kFALSE),
    fZenithMin(0),
    fZenithMax(0),
    fAzimuthMin(0),
    fAzimuthMax(0),
    fPRange(0),
    fPResolution(1),
    fAp(0),
    fMomentumDist(0),
    fUnfoldedMomentumDist(0),
    fZenithDist(0),
    fPDist(0),
    fNParticles(0)
{
  //
  // Standard ctor.
  //
  fName = "ACORDE";
  fTitle = "Cosmic Muons generator";

  // Set the origin above the vertex, on the surface.
  fOrigin[0] = 0.;
  fOrigin[1] = AliACORDEConstants::Instance()->Depth(); // At the surface by default.
  fOrigin[2] = 0.;
}

//_____________________________________________________________________________
AliGenACORDE::~AliGenACORDE()
{
  //
  // Default dtor.
  //
  if ( fPDist ) {fPDist->Delete(); delete fPDist; fPDist = 0;}
  if ( fUnfoldedMomentumDist ) delete fUnfoldedMomentumDist;
  if ( fMomentumDist ) delete fMomentumDist;
  if ( fAp )           delete fAp;
  if ( fCRModeName )   delete fCRModeName;
}

//_____________________________________________________________________________
void AliGenACORDE::Generate()
{
  //
  // Generate on one trigger
  // Call the respective method inside the loop for the number
  // of tracks per trigger.

  for (Int_t i = 0; i < fNParticles; i++ ) {

    if ( fCRMode == kMuonBundle ) {
      this->GenerateOneMuonBundle();

    } else if ( fCRMode == kSingleMuons ) {
      this->GenerateOneSingleMuon(kTRUE);

    } else {
      // Generate only single muons following the parametrizations
      // for atmospheric muons.
      this->GenerateOneSingleMuon();

    }

  }
}

//_____________________________________________________________________________
void AliGenACORDE::Init()
{
  //
  // Initialize some internal methods.
  //


   printf("**************************************************************\n");
   printf("<<< *** Starting the AliGenACORDE cosmic generator ******** >>>\n");
   printf("<<< *** No. of muons generated at the surface of P2: %d,  * >>>\n",fNParticles);
   printf("**************************************************************\n");

  // Determine some specific data members.
  fPRange = TMath::Abs(fPMax-fPMin);

  if ( fCRMode == kSingleMuons ) {
    fCRModeName = new TString("Single Muons");
    // Initialisation, check consistency of selected ranges
    if(TestBit(kPtRange)&&TestBit(kMomentumRange)) 
      Fatal("Init","You should not set the momentum range and the pt range!");
    
    if((!TestBit(kPtRange))&&(!TestBit(kMomentumRange))) 
      Fatal("Init","You should set either the momentum or the pt range!");
    
  } else if ( fCRMode == kMuonBundle ) {
    fCRModeName = new TString("Muon Bundles");

  } else if ( fCRMode == kMuonFlux ) {
    fCRModeName = new TString("Muon Fluxes");
    // Initialize the ditribution functions.
    this->InitMomentumGeneration();
    this->InitZenithalAngleGeneration();
    
  } else {
    Fatal("Generate", "Generation Mode unknown!\n");

  }

}

//____________________________________________________________________________
void AliGenACORDE::GenerateOneSingleMuon(Bool_t withFlatMomentum)
{
  //
  // Generate One Single Muon
  // This method will generate only one muon.
  // The momentum will be randomly flat distributed if
  // the paremeter "withFlatMomentum" is set to kTRUE,
  // otherwise the momentum will generate acordingly the parametrization
  // given by 
  // and adpted from Bruno Alessandro's implementation with the
  // CERNLIB to AliRoot.
  // The "withFlatMomentum" parameter also will be used to generate
  // the muons with a flat Zenithal angle distribution.
  // Do the smearing here, so that means per track.

  Float_t polar[3]= {0,0,0}; // Polarization parameters
  Float_t origin[3];
  Int_t nt;
  Float_t p[3];
  Float_t pmom, pt;
  Float_t random[6];

  // Take the azimuth random.
  Rndm(random, 2);
  Float_t azimuth = fAzimuthMin + (fAzimuthMax-fAzimuthMin)*random[0];
  Float_t zenith = fZenithMin + (fZenithMax - fZenithMin)*random[1];

  if ( withFlatMomentum ) {
    Rndm(random,3);
    if(TestBit(kMomentumRange)) {
      pmom = -( fPMin + random[0]*(fPMax - fPMin) ); // always downwards.
      pt = pmom*TMath::Sin(zenith*kDegrad);
    } else {
      pt = -( fPtMin + random[1]*(fPtMax - fPtMin)); // always downwards.
      pmom = pt/TMath::Sin(zenith*kDegrad);
    }

  } else {
    if ( fMomentumDist ) {
      pmom = -this->GetMomentum(); // Always downwards.
    } else {
      pmom = -fPMin;
    }
    zenith = this->GetZenithAngle(pmom);  // In degrees
    pt = pmom*TMath::Sin(zenith*kDegrad);
  }

  p[0] = pt*TMath::Sin(azimuth*kDegrad);
  p[1] = pmom*TMath::Cos(zenith*kDegrad);
  p[2] = pt*TMath::Cos(azimuth*kDegrad);

  // Finaly the origin, with the smearing
  Rndm(random,6);
  origin[0] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
    TMath::Sin(azimuth*kDegrad)
    + fOsigma[0]* TMath::Cos(2*random[0]*TMath::Pi())*TMath::Sqrt(-2*TMath::Log(random[1]));

  origin[1] = AliACORDEConstants::Instance()->Depth();

  origin[2] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
    TMath::Cos(azimuth*kDegrad)
    + fOsigma[2]* TMath::Cos(2*random[2]*TMath::Pi())*TMath::Sqrt(-2*TMath::Log(random[3]));

  // Put the track on the stack.
  PushTrack(fTrackIt,-1,fIpart,p,origin,polar,0,kPPrimary,nt);

}

//____________________________________________________________________________
void AliGenACORDE::GenerateOneMuonBundle()
{
  //
  // Generate One Muon Bundle method
  // This method will generate a bunch of muons following the
  // procedure of the AliGenScan class.
  // These muons will be generated with flat momentum.

  Float_t polar[3]= {0,0,0}; // Polarization parameters
  Float_t origin[3];
  Float_t p[3];
  Int_t nt;
  Float_t pmom;
  Float_t random[6];

  Rndm(random, 3);
  Float_t zenith = fZenithMin + (fZenithMax - fZenithMin)*random[1];
  Float_t azimuth = fAzimuthMin + (fAzimuthMax-fAzimuthMin)*random[2];
  //Float_t zenith = 10;
  //Float_t azimuth = 30;

  // Generate the kinematics a la AliGenScan (Andreas Morchs)
  Float_t dx, dz;
  if ( fNx > 0 ) {
    dx = fXwidth/fNx;
  } else {
    dx = 1e10;
    //dx = 100.;
  }

  if ( fNz > 0 ) {
    dz = fZwidth/fNz;
  } else {
    dz = 1e10;
    //dz = 100.;
  }

  origin[0] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
              TMath::Sin(azimuth*kDegrad);
  //origin[0] = 0.;
  origin[1] = AliACORDEConstants::Instance()->Depth();
  //origin[1] = 900;
  origin[2] = AliACORDEConstants::Instance()->Depth()*TMath::Tan(zenith*kDegrad)*
              TMath::Cos(azimuth*kDegrad);
    //origin[2] = 0.;

  for (Int_t ix = 0; ix < fNx; ix++ ) {
    for (Int_t iz = 0; iz < fNz; iz++ ) {
      Rndm(random,6);
      origin[0]+=ix*dx+2*(random[1]-0.5)*fOsigma[0];
      origin[2]+=iz*dz+2*(random[2]-0.5)*fOsigma[2];
      if ( random[4] < 0.5 ) {
        origin[0] = -1*origin[0];
      }
      if ( random[5] < 0.5 ) {
        origin[2] = -1*origin[2];
      }

      pmom = -(fPMin + random[3] *(fPMax - fPMax) ); // Always downwards
      p[0] = TMath::Sin(zenith*kDegrad)*TMath::Sin(azimuth*kDegrad)*pmom;
      p[1] = TMath::Cos(zenith*kDegrad)*pmom;
      p[2] = TMath::Sin(zenith*kDegrad)*TMath::Cos(azimuth*kDegrad)*pmom;

      PushTrack(fTrackIt, -1, fIpart, p, origin, polar, 0, kPPrimary, nt);
    }

  }

}

//____________________________________________________________________________
void AliGenACORDE::SetGridRange(Int_t nx,Float_t xwidth, Int_t nz, Float_t zwidth)
{
  //
  // Define the grid
  // This data shuold be used for Muon bundles generation.
  //
  fXwidth=xwidth;
  fNx=nx;
  fZwidth=zwidth;
  fNz=nz;

  // Print a message  about the use, if the Mode has not been set, or
  // it has to a different Mode.
  if ( fCRMode != kMuonBundle ) {
    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");
    fMuonGrid = kTRUE;
  }
}

//____________________________________________________________________________
void AliGenACORDE::InitApWeightFactors()
{
  //
  // This factors will be  to correct the zenithal angle distribution
  // acording the momentum

  //
  // Fill the array for the flux zenith angle dependence.
  // at the index 0 of fAp[] will be the "factor" if we have a muon
  // of 0 GeV.
  Float_t a[6] = {-1.61, -1.50, -1.28, -0.94, -0.61, -0.22};
  Float_t p[6] = { 0., 10., 30., 100., 300., 1000.};

  // Get the information from the momentum
  Int_t pEnd  = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
  // Initialize the Array of floats to hold the a(p) factors.
  fAp = new TArrayF(pEnd);
  
  Int_t index = 0;

  for (Int_t i = 0; i < pEnd; i++ ) {
    Float_t currentP = ((Float_t)i)*fPResolution;
    if ( currentP < p[1] )                          index = 0;
    else if ( currentP >= p[1] && currentP < p[2] ) index = 1;
    else if ( currentP >= p[2] && currentP < p[3] ) index = 2;
    else if ( currentP >= p[3] && currentP < p[4] ) index = 3;
    else if ( currentP >= p[4] )                    index = 4;

    Float_t ap = (currentP -p[index])*(a[index+1] - a[index])/
                 (p[index+1] - p[index]) + a[index];
    fAp->AddAt(ap, i);
  }

}

//___________________________________________________________________________
void AliGenACORDE::InitMomentumGeneration()
{
  //
  // Initialize a funtion to generate the momentum randomly
  // acording this function.
  //

  // Check if we nned to initialize the function
  if ( fPMin != fPMax ) {

    // Check also if the function have been defined yet.
    if ( !fMomentumDist ) {

      // If not, use this function
      const char* y      = "log10(x)";
      
      const char* h1Coef = "[0]*( %s*%s*%s/2 - (5*%s*%s/2) + 3*%s )";
      const char* h2Coef = "[1]*( (-2*%s*%s*%s/3) + (3*%s*%s) - 10*%s/3 + 1 )";
      const char* h3Coef = "[2]*( %s*%s*%s/6 - %s*%s/2 + %s/3 )";
      const char* s2Coef = "[3]*( %s*%s*%s/3 - 2*%s*%s + 11*%s/3 - 2 )";
      
      const char* h = "%s + %s + %s + %s";
      const char* flux = "pow(10., %s)";
      const char* normalizedFlux = "0.86*x*x*x*pow(10., %s)";
      const char* paramNames[4] = {"H1", "H2", "H3", "S1"};
      
      char buffer1[1024];
      char buffer2[1024];
      char buffer3[1024];
      char buffer4[1024];
      char buffer5[1024];
      char buffer6[1024];
      char buffer7[1024];

      sprintf(buffer1, h1Coef, y, y, y, y, y, y);
      sprintf(buffer2, h2Coef, y, y, y, y, y, y);
      sprintf(buffer3, h3Coef, y, y, y, y, y, y);
      sprintf(buffer4, s2Coef, y, y, y, y, y, y);
      
      sprintf(buffer5, h, buffer1, buffer2, buffer3, buffer4);
      
      sprintf(buffer6, flux, buffer5);
      
      fMomentumDist = new TF1("fMomentumDist", buffer6, fPMin, fPMax);
      sprintf(buffer7, normalizedFlux, buffer5);
      fUnfoldedMomentumDist = new TF1("fUnfoldedMomentumDist", buffer7, fPMin, fPMax);
      for (Int_t i = 0; i < 4; i++ ) {
	fMomentumDist->SetParName(i, paramNames[i]);
	fUnfoldedMomentumDist->SetParName(i, paramNames[i]);
      }
      
      fMomentumDist->SetParameter(0, 0.133);
      fMomentumDist->SetParameter(1, -2.521);
      fMomentumDist->SetParameter(2, -5.78);
      fMomentumDist->SetParameter(3, -2.11);

      fUnfoldedMomentumDist->SetParameter(0, 0.133);
      fUnfoldedMomentumDist->SetParameter(1, -2.521);
      fUnfoldedMomentumDist->SetParameter(2, -5.78);
      fUnfoldedMomentumDist->SetParameter(3, -2.11);
      
    }

  }

}

//____________________________________________________________________________
void AliGenACORDE::InitZenithalAngleGeneration()
{
  //
  // Initalize a distribution function for the zenith angle.
  // This angle will be obtained randomly acording this function.
  // The generated angles  will been in degrees.

  // Check if we need to create the function.
  if ( fZenithMin != fZenithMax ) {

    // Check also if another function have been defined.
    if ( !fZenithDist ) {
      
      // initialize the momentum dependent coefficients, a(p) 
      this->InitApWeightFactors();

      Int_t pEnd  = TMath::Abs(TMath::Nint(fPMax/fPResolution)) + 1;
      char name[26];
      char title[52];
      fPDist = new TClonesArray("TH1F", pEnd);
      TClonesArray &mom = *fPDist;
      TH1F* zenith = 0;
      Float_t weight = 0;
      for ( Int_t i = 0; i < pEnd; i++ ) {
	// Fill the distribution
	sprintf(name, "zenith%d", i+1);
	sprintf(title, "Zenith distribution, p=%f", fPMin+(Float_t)i);
	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));

	// Make a loop for the angle and fill the histogram for the weight
	Int_t steps = 1000;
	Float_t value = 0;
	for (Int_t j = 0; j < steps; j++ ) {
	  value = TMath::Cos(fZenithMin*TMath::Pi()/180) + (Float_t)j * ( TMath::Cos(fZenithMax*TMath::Pi()/180) - TMath::Cos(fZenithMin*TMath::Pi()/180))/1000;
	  weight = 1 + fAp->At(i)*(1 - value);
	  zenith->Fill(value, weight);
	}

      }

    } 

  }

}

//____________________________________________________________________________
Float_t AliGenACORDE::GetZenithAngle(Float_t mom) const
{

  Float_t zenith = 0.;
  // Check if you need to generate a constant zenith angle.
  if ( !fZenithDist ) {
    // Check if you have defined an array of momentum functions
    if ( fPDist ) {
      Int_t pIndex = TMath::Abs(TMath::Nint(mom));
      TH1F* cosZenithAngle = (TH1F*)fPDist->UncheckedAt(pIndex);
      Float_t tmpzenith = TMath::ACos(cosZenithAngle->GetRandom());
      // Correct the value
      zenith = kRaddeg*tmpzenith;
      return zenith;
    } else {

      if ( fCRMode != kMuonFlux ) {
	// If you aren't generating muons obeying any ditribution
	// only generate a flat zenith angle, acording the input settings
	Float_t random[2];
	Rndm(random, 2);
	zenith = fZenithMin + (fZenithMax - fZenithMin)*random[0];

      } else {
	// Even if you are generating muons acording some distribution,
	// but you don't want to ...
	zenith = fZenithMin;

      }

    }
  } else {
    zenith = fZenithDist->GetRandom();
  }

  return zenith;
}

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