EVGEN/AliGenMUONLMR.cxx

   1 //#include <TClonesArray.h>
   2
   3 #include <TDatabasePDG.h>
   4 #include <TFile.h>
   5 #include "AliConst.h"
   6 #include "AliGenMUONLMR.h"
   7 #include "AliMC.h"
   8 #include "AliRun.h"
   9 #include "AliLog.h"
  10 #include "AliGenEventHeader.h"
  11
  12 ClassImp(AliGenMUONLMR)
  13
  14   AliGenMUONLMR::AliGenMUONLMR () : AliGenMC(),
  15                                     fNMuMin(2),
  16                                     fCMSEnergy(kNCMSEnergies),
  17                                     fGenSingleProc(-1),
  18                                     fYCM(0),
  19                                     fCosTheta(0x0),
  20                                     fRhoLineShape(0x0),
  21                                     fHMultMu(0x0),
  22                                     fHNProc(0x0) {
  23     //
  24     // default constructor
  25     //
  26     for (int i=0; i<fgkNpart; i++) {
  27       fPDG[i] = 0;
  28       fScaleMult[i] = 1.;
  29       fPt[i] = NULL;
  30       fY[i] = NULL;
  31       fMult[i] = NULL;
  32       fParticle[i] = NULL;
  33     }
  34     for (int i=0; i<2; i++) {
  35       fMu[i]     = NULL;
  36       fDecay[i] = NULL;
  37     }
  38
  39     for (int i=0; i<3; i++) {
  40       fDalitz[i] = NULL;
  41     }
  42
  43   }
  44
  45 //-----------------------------------------------------------
  46
  47 void AliGenMUONLMR::SetCMSEnergy(CMSEnergies energy){
  48   fCMSEnergy = energy;
  49   // initialize pt and y distributions according to a fit to
  50   // Pythia simulation at sqrt(s) = 7 TeV
  51   for (Int_t ipart=0; ipart < fgkNpart; ipart++) fScaleMult[ipart] = 1;
  52   fScaleMult[kPionLMR] = 0; // set pion multiplicity to zero
  53   fScaleMult[kKaonLMR] = 0; // set kaon multiplicity to zero
  54   const char* fdname[2] = {"fDecPion","fDecKaon"};
  55   Double_t ctau[2] = {7.8045, 3.712};
  56   Int_t pdg[7] = {211, 321, 221, 113, 223, 333, 331};
  57   const char* fptname[7] = {"fPtPion","fPtKaon","fPtEta","fPtRho","fPtOmega","fPtPhi","fPtEtaPrime"};
  58   const char* fyname[7] = {"fYPion","fYKaon","fYEta","fYRho","fYOmega","fYPhi","fYEtaPrime"};
  59   const char* fnname[7] = {"fMultPion","fMultKaon","fMultEta","fMultRho","fMultOmega","fMultPhi","fMultEtaPrime"};
  60   Double_t ptparam[7][9];
  61   Double_t yparam[7][9];
  62   Double_t nparam[7][9];
  63
  64   // parameters for 8 TeV generation
  65   if (fCMSEnergy==kCMS8000GeV) {
  66     AliInfo ("Using pp parameterization at 8 TeV\n");
  67
  68   // Parameters of transverse momentum spectra
  69   Double_t ptparam8000[7][9] = {{1,0.427,2.52,0,0,0,0,0,0},   // pions from 7 TeV section of code
  70                                 {1,0.58,2.57,0,0,0,0,0,0},    // kaons from 7 TeV section of code
  71                                 {1,0.657,2.685,0,0,0,0,0,0},  // eta from PYTHIA6.4 ATLAS-CSC at 8 TeV
  72                                 {1,1.44,3.16,0,0,0,0,0,0},    // rho+omega from 7 TeV section of code
  73                                 {1,1.44,3.16,0,0,0,0,0,0},    // rho+omega from 7 TeV section of code
  74                                 {1,1.16,2.74,0,0,0,0,0,0},    // phi from 7 TeV section of code
  75                                 {1,0.755,2.578,0,0,0,0,0,0}}; // etaPrime from PYTHIA 6.4 ATLAS-CSC at 8 TeV
  76
  77   // Parameters of rapidity spectra
  78   Double_t yparam8000[7][9] = {{1,0.8251,3.657,0,0,0,0,0,0},  // pions from 7 TeV section of code
  79                                {1,1.83,2.698,0,0,0,0,0,0},    // kaons from 7 TeV section of code
  80                                {1,0.0509,3.96,0,0,0,0,0,0},   // eta from PYTHIA6.4 ATLAS-CSC at 8 TeV
  81                                {1,0.0489,3.961,0,0,0,0,0,0},  // rho from PYTHIA6.4 ATLAS-CSC at 8 TeV
  82                                {1,0.0650,3.966,0,0,0,0,0,0},  // omega from PYTHIA 6.4 ATLAS-CSC at 8 TeV
  83                                {1,1.279,2.745,0,0,0,0,0,0},   // phi from from PYTHIA6.4 ATLAS-CSC at 8 TeV
  84                                {1,0.1627,3.883,0,0,0,0,0,0}}; // eta prime from PYTHIA6.4 ATLAS-CSC at 8 TeV
  85
  86   // Parameters of multiplicity spectra
  87   Double_t nparam8000[7][9] = {{353.582, 6.76263, 1.66979, 998.445, 9.73281, 12.6704, 175.187, 29.08, 40.2531}, //pions from 7 TeV section of code
  88                                {1.e4,    0.2841, 0,0,0,0,0,0,0},     // kaons from 7 TeV section of code
  89                                {2.279e4, 0.2622, 0,0,0,0,0,0,0},  // eta from PYTHIA6.4 ATLAS-CSC at 8 TeV
  90                                {1.564e4, 0.1713, 0,0,0,0,0,0,0},  // rho from PYTHIA6.4 ATLAS-CSC at 8 TeV
  91                                {1.662e4, 0.183,  0,0,0,0,0,0,0},   // omega from PYTHIA6.4 ATLAS-CSC at 8 TeV
  92                                {6.723e4, 1.121,  0,0,0,0,0,0,0},  // phi from PYTHIA6.4 ATLAS-CSC at 8 TeV
  93                                {5.005e4, 0.6971, 0,0,0,0,0,0,0}}; // eta prime from PYTHIA6.4 ATLAS-CSC at 8 TeV
  94
  95     for (Int_t i=0; i<fgkNpart; i++) {
  96       for (Int_t j=0; j<9; j++) {
  97         ptparam[i][j] = ptparam8000[i][j];
  98         yparam[i][j] = yparam8000[i][j];
  99         nparam[i][j] = nparam8000[i][j];
 100       }
 101     }
 102   }
 103
 104   // parameters for 7 TeV generation
 105   else if (fCMSEnergy==kCMS7000GeV) {
 106     AliInfo ("Using pp parameterization at 7 TeV\n");
 107     Double_t ptparam7000[7][9] = {{1,0.427,2.52,0,0,0,0,0,0}, // pions from Pythia
 108                                   {1,0.58,2.57,0,0,0,0,0,0},  // kaons from Pythia
 109                                   {1,0.641,2.62,0,0,0,0,0,0}, // eta from Pythia
 110                                   {1,1.44,3.16,0,0,0,0,0,0},  // rho+omega from ALICE muon
 111                                   {1,1.44,3.16,0,0,0,0,0,0},  // rho+omega from ALICE muon
 112                                   {1,1.16,2.74,0,0,0,0,0,0},  // phi from ALICE muon
 113                                   {1,0.72,2.5,0,0,0,0,0,0}};  // etaPrime from Pythia
 114
 115     Double_t yparam7000[7][9] = {{1,0.8251,3.657,0,0,0,0,0,0}, // pions from pythia
 116                                  {1,1.83,2.698,0,0,0,0,0,0},   // kaons from pythia
 117                                  {1,1.169,3.282,0,0,0,0,0,0},  // eta from pythia
 118                                  {1,1.234,3.264,0,0,0,0,0,0},  // rho from pythia
 119                                  {1,1.311,3.223,0,0,0,0,0,0},  // omega from pythia
 120                                  {1,2.388,2.129,0,0,0,0,0,0},  // phi from pythia
 121                                  {1,1.13,3.3,0,0,0,0,0,0}};    // eta prime from pythia
 122
 123     // multiplicity parameters from pythia
 124     Double_t nparam7000[7][9] = {{353.582, 6.76263, 1.66979, 998.445, 9.73281, 12.6704, 175.187, 29.08, 40.2531},
 125                                  {1.e4,  0.2841, 0,0,0,0,0,0,0},
 126                                  {1.e4,  0.2647, 0,0,0,0,0,0,0},
 127                                  {7055,  0.1786, 0,0,0,0,0,0,0},
 128                                  {7500,  0.1896, 0,0,0,0,0,0,0},
 129                                  {5.e4,  1.167,  0,0,0,0,0,0,0},
 130                                  {2.9e4, 0.714,  0,0,0,0,0,0,0}};
 131
 132     for (Int_t i=0; i<fgkNpart; i++) {
 133       for (Int_t j=0; j<9; j++) {
 134         ptparam[i][j] = ptparam7000[i][j];
 135         yparam[i][j] = yparam7000[i][j];
 136         nparam[i][j] = nparam7000[i][j];
 137       }
 138     }
 139   }
 140   else if (fCMSEnergy==kCMS5020GeVpPb || fCMSEnergy==kCMS5020GeVPbp) {
 141     AliInfo ("Using pPb parameterization at 5.02 TeV\n");
 142     Double_t ptparam5020[7][9] = {{1,0.427,2.52,0,0,0,0,0,0}, // pions from Pythia at 7 TeV
 143                                   {1,0.58,2.57,0,0,0,0,0,0},  // kaons from Pythia at 7 TeV
 144                                   {1,0.665,2.796,0,0,0,0,0,0}, // eta from Pythia at 5.02 TeV
 145                                   {1,1.66,3.12,0,0,0,0,0,0},  // rho+omega from ALICE muon
 146                                   {1,1.66,3.12,0,0,0,0,0,0},  // rho+omega from ALICE muon
 147                                   {1,2.03,3.13,0,0,0,0,0,0},  // phi from ALICE muon
 148                                   {1,0.767,2.713,0,0,0,0,0,0}};  // etaPrime from Pythia at 5.02 TeV
 149
 150     Double_t yparam5020[7][9] = {{1,0.8251,3.657,0,0,0,0,0,0}, // pions from pythia at 7 TeV
 151                                  {1,1.83,2.698,0,0,0,0,0,0},   // kaons from pythia at 7 TeV
 152                                  {1,1.169,3.282,0,0,0,0,0,0},  // eta from pythia at 7 TeV
 153                                  {1,1.234,3.264,0,0,0,0,0,0},  // rho from pythia at 7 TeV
 154                                  {1,1.311,3.223,0,0,0,0,0,0},  // omega from pythia at 7 TeV
 155                                  {1,2.388,2.129,0,0,0,0,0,0},  // phi from pythia at 7 TeV
 156                                  {1,1.13,3.3,0,0,0,0,0,0}};    // eta prime from pythia at 7 TeV
 157
 158     // multiplicity parameters from pythia at 7 TeV
 159     Double_t nparam5020[7][9] = {{353.582, 6.76263, 1.66979, 998.445, 9.73281, 12.6704, 175.187, 29.08, 40.2531},
 160                                  {1.e4,  0.2841, 0,0,0,0,0,0,0},
 161                                  {1.e4,  0.2647, 0,0,0,0,0,0,0},
 162                                  {7055,  0.1786, 0,0,0,0,0,0,0},
 163                                  {7500,  0.1896, 0,0,0,0,0,0,0},
 164                                  {5.e4,  1.167,  0,0,0,0,0,0,0},
 165                                  {2.9e4, 0.714,  0,0,0,0,0,0,0}};
 166
 167     for (Int_t i=0; i<fgkNpart; i++) {
 168       for (Int_t j=0; j<9; j++) {
 169         ptparam[i][j] = ptparam5020[i][j];
 170         yparam[i][j] = yparam5020[i][j];
 171         nparam[i][j] = nparam5020[i][j];
 172       }
 173     }
 174     if (fCMSEnergy==kCMS5020GeVpPb) fYCM = -0.4654;
 175     else fYCM = 0.4654;
 176   }
 177   else if (fCMSEnergy==kCMS2760GeV){
 178     // parameters for 2.76 generation
 179     // pt params has been determined as <pt>ALICE_2.76 = <pt>ALICE_7 * <pt>PYTHIA_2.76 / <pt>PYTHIA_7
 180     AliInfo ("Using pp parameterization at 2.76 TeV\n");
 181     Double_t yparam2760[7][9] = {{1,0.8251,3.657,0,0,0,0,0,0},      // pions from pythia
 182                                  {1,1.83,2.698,0,0,0,0,0,0},        // kaons from pythia
 183                                  {1,0.011,3.474,0,0,0,0,0,0},       // eta from pythia
 184                                  {1,-0.01,3.409,0,0,0,0,0,0},       // rho from pythia
 185                                  {1,-0.037,3.294,0,0,0,0,0,0},      // omega from pythia
 186                                  {1,-0.016,2.717,0,0,0,0,0,0},      // phi from pythia
 187                                  {1,-0.010,3.312,0,0,0,0,0,0}};     // eta prime from pythia
 188
 189     Double_t ptparam2760[7][9] = {{1,0.1665,8.878,0,0,0,0,0,0},  // pions from Pythia
 190                                   {1,0.1657,8.591,0,0,0,0,0,0},  // kaons from Pythia
 191                                   {1,0.641,2.62,0,0,0,0,0,0},    // eta from ALICE 7 TeV
 192                                   {1,1.3551,3.16,0,0,0,0,0,0},   // rho with <pt> scaled
 193                                   {1,1.3551,3.16,0,0,0,0,0,0},   // omega with <pt> scaled
 194                                   {1,1.0811,2.74,0,0,0,0,0,0},   // phi with <pt> scaled
 195                                   {1,0.72,2.5,0,0,0,0,0,0}};     // etaPrime from ALICE 7 TeV
 196
 197     Double_t nparam2760[7][9] = {{9752,-2.693,3.023,9.5e9,-84.68,16.75,-14.06,635.3,-423.2}, // pions
 198                                  {1.e5, 1.538,  0,0,0,0,0,0,0},                              // kaons
 199                                  {1.e4, 0.351,  0,0,0,0,0,0,0},                              // eta
 200                                  {1.e4, 0.2471, 0,0,0,0,0,0,0},                              // rho
 201                                  {1.e4, 0.2583, 0,0,0,0,0,0,0},                              // omega
 202                                  {1.e5, 1.393,  0,0,0,0,0,0,0},                              // phi
 203                                  {1.e4, 0.9005, 0,0,0,0,0,0,0}};                             // etaPrime
 204
 205     for (Int_t i=0; i<fgkNpart; i++) {
 206       for (Int_t j=0; j<9; j++) {
 207         ptparam[i][j] = ptparam2760[i][j];
 208         yparam[i][j] = yparam2760[i][j];
 209         nparam[i][j] = nparam2760[i][j];
 210       }
 211     }
 212   }
 213   else AliFatal("Energy not correctly defined");
 214
 215   for (Int_t i=0; i<fgkNpart; i++) {
 216     fPDG[i] = pdg[i];
 217     if (i!=0) {
 218       fMult[i] = new TF1(fnname[i],"[0]*exp(-[1]*x)",0,30);
 219       fMult[i]->SetParameters(nparam[i][0],nparam[i][1]);
 220     }
 221     else {
 222       fMult[i] = new TF1(fnname[i],"gaus(0)+gaus(3)+gaus(6)",0,150);
 223       for (Int_t j=0; j<9; j++) fMult[i]->SetParameter(j,nparam[i][j]);
 224     }
 225
 226     fPt[i] = new TF1(fptname[i],AliGenMUONLMR::PtDistr,0,20,3);
 227     fPt[i]->SetParameters(ptparam[i][0], ptparam[i][1], ptparam[i][2]);
 228     fY[i] = new TF1(fyname[i],AliGenMUONLMR::YDistr,-10,10,4);
 229     fY[i]->SetParameters(yparam[i][0], yparam[i][1], yparam[i][2],fYCM);
 230   }
 231
 232   for(Int_t i = 0; i<2; i++){
 233     fDecay[i] = new TF1(fdname[i],"exp(-x/[0])",0,150);
 234     fDecay[i]->SetParameter(0,ctau[i]);
 235   }
 236
 237   for (Int_t ipart = 0; ipart < fgkNpart; ipart++) {
 238     fParticle[ipart] = new TParticle();
 239     fParticle[ipart]->SetPdgCode(fPDG[ipart]);
 240   }
 241
 242   TDatabasePDG *pdgdb = TDatabasePDG::Instance();
 243   Double_t mumass = pdgdb->GetParticle(13)->Mass();
 244   fMu[0] = new TParticle();
 245   fMu[0]->SetPdgCode(-13);
 246   fMu[0]->SetCalcMass(mumass);
 247   fMu[1] = new TParticle();
 248   fMu[1]->SetPdgCode(13);
 249   fMu[1]->SetCalcMass(mumass);
 250
 251   // function for polarized theta distributions
 252   fCosTheta = new TF1 ("fCosTheta","1+[0]*x*x",-1,1);
 253   fCosTheta->SetParameter(0,1);
 254
 255   // Dalitz decays
 256   Int_t nbins = 1000;
 257   Double_t xmin = 0, xmax = 2;
 258   fDalitz[0] = new TH1F("hDalitzEta","",nbins,xmin,xmax);
 259   fDalitz[1] = new TH1F("hDalitzOmega","",nbins,xmin,xmax);
 260   fDalitz[2] = new TH1F("hDalitzEtaPrime","",nbins,xmin,xmax);
 261
 262   Double_t meta   = pdgdb->GetParticle("eta")->Mass();
 263   Double_t momega = pdgdb->GetParticle("omega")->Mass();
 264   Double_t metaPrime = pdgdb->GetParticle("eta'")->Mass();
 265   Double_t mpi0   = pdgdb->GetParticle("pi0")->Mass();
 266   Double_t md3 = 0, mres = 0;
 267
 268   for (Int_t index = 0; index < 3; index++) {
 269     if (index == 0) {
 270       mres = meta;
 271       md3 = 0;
 272     }
 273     else if (index == 1) {
 274       mres = momega;
 275       md3 = mpi0;
 276     }
 277     else if (index == 2) {
 278       mres = metaPrime;
 279       md3 = 0;
 280     }
 281     Double_t delta   = md3 * md3 / (mres * mres);
 282     Double_t epsilon = mumass * mumass / (mres * mres);
 283     nbins = fDalitz[index]->GetNbinsX();
 284     xmin = fDalitz[index]->GetXaxis()->GetXmin();
 285     Double_t deltax =  fDalitz[index]->GetBinWidth(1);
 286     Double_t xd = xmin - deltax/2.;
 287     for (Int_t ibin = 0; ibin< nbins; ibin++) {
 288       Double_t dalval = 0;
 289       xd += deltax;
 290       if (xd > 4. *epsilon) {
 291         Double_t bracket = TMath::Power(1. + xd/(1. - delta),2)
 292           - 4. * xd / ((1. - delta) * (1. - delta));
 293         if (bracket > 0) {
 294           dalval = TMath::Power(bracket,1.5) /xd *
 295             TMath::Sqrt(1 - 4 * epsilon / xd) * (1 + 2 * epsilon / xd) *
 296             FormFactor(xd * mres * mres, index);
 297           fDalitz[index]->Fill(xd,dalval);
 298         }
 299       }
 300     }
 301   }
 302
 303   fRhoLineShape = new TF1("fRhoLineShape",RhoLineShapeNew,0,2,2);
 304   fHMultMu = new TH1D("fHMultMu","Muon multiplicity",20,-0.5,19.5);
 305   fHNProc = new TH1D("fHNProc","Number of gen. evts. per process in 4 pi",9,-0.5,8.5);
 306
 307 }
 308 //-----------------------------------------------------------
 309
 310 AliGenMUONLMR::AliGenMUONLMR (AliGenMUONLMR &gen) : AliGenMC(),
 311                                                     fNMuMin(gen.fNMuMin),
 312                                                     fCMSEnergy(gen.fCMSEnergy),
 313                                                     fGenSingleProc(gen.fGenSingleProc),
 314                                                     fYCM(gen.fYCM),
 315                                                     fCosTheta(gen.fCosTheta),
 316                                                     fRhoLineShape(gen.fRhoLineShape),
 317                                                     fHMultMu(gen.fHMultMu),
 318                                                     fHNProc(gen.fHNProc) {
 319   for (Int_t i=0; i < fgkNpart; i++) {
 320     fPDG[i] = gen.fPDG[i];
 321     fScaleMult[i] = gen.fScaleMult[i];
 322     fPt[i] = (TF1*) gen.fPt[i]->Clone();
 323     fY[i] = (TF1*) gen.fY[i]->Clone();
 324     fMult[i] = (TF1*) gen.fMult[i]->Clone();
 325     fParticle[i] = (TParticle*) gen.fParticle[i]->Clone();
 326   }
 327
 328   for(Int_t i = 0; i<2; i++) fDecay[i] = (TF1*) gen.fDecay[i]->Clone();
 329   for(Int_t i = 0; i<3; i++) fDalitz[i] = (TH1F*) gen.fDalitz[i]->Clone();
 330   for(Int_t i = 0; i<2; i++) fMu[i] = (TParticle*) gen.fMu[i]->Clone();
 331 }
 332
 333 //-----------------------------------------------------------
 334
 335 AliGenMUONLMR& AliGenMUONLMR::operator=(const AliGenMUONLMR &gen) {
 336   // Assignment operator
 337   if (this!=&gen) {
 338     fNMuMin = gen.fNMuMin;
 339     fCMSEnergy = gen.fCMSEnergy;
 340     fGenSingleProc = gen.fGenSingleProc;
 341     fYCM = gen.fYCM;
 342     fCosTheta = (TF1*) gen.fCosTheta->Clone();
 343     fRhoLineShape = (TF1*) gen.fRhoLineShape->Clone();
 344     fHMultMu = (TH1D*) gen.fHMultMu->Clone();
 345     fHNProc = (TH1D*) gen.fHNProc->Clone();
 346
 347     for (Int_t i=0; i < fgkNpart; i++) {
 348       fPDG[i] = gen.fPDG[i];
 349       fScaleMult[i] = gen.fScaleMult[i];
 350       fPt[i] = (TF1*) gen.fPt[i]->Clone();
 351       fY[i] = (TF1*) gen.fY[i]->Clone();
 352       fMult[i] = (TF1*) gen.fMult[i]->Clone();
 353       fParticle[i] = (TParticle*) gen.fParticle[i]->Clone();
 354     }
 355
 356     for(Int_t i = 0; i<2; i++) fDecay[i] = (TF1*) gen.fDecay[i]->Clone();
 357     for(Int_t i = 0; i<3; i++) fDalitz[i] = (TH1F*) gen.fDalitz[i]->Clone();
 358     for(Int_t i = 0; i<2; i++) fMu[i] = (TParticle*) gen.fMu[i]->Clone();
 359   }
 360   return *this;
 361 }
 362
 363
 364 //-----------------------------------------------------------
 365
 366 AliGenMUONLMR::~AliGenMUONLMR()
 367 {
 368   // Default destructor
 369   for (Int_t i=0; i<7; i++) {
 370     delete fPt[i];
 371     delete fY[i];
 372     delete fMult[i];
 373     delete fParticle[i];
 374   }
 375
 376   for (Int_t i=0; i<2; i++) {
 377     delete fDecay[i];
 378     delete fMu[i];
 379   }
 380
 381   for (Int_t i=0; i<3; i++) delete fDalitz[i];
 382
 383   delete fCosTheta; fCosTheta = 0;
 384   delete fRhoLineShape; fRhoLineShape = 0;
 385   delete fHMultMu; fHMultMu = 0;
 386   delete fHNProc;  fHNProc = 0;
 387 }
 388
 389 //-----------------------------------------------------------
 390
 391 void AliGenMUONLMR::FinishRun(){
 392   // save some histograms to an output file
 393   Int_t nbins = fHNProc->GetNbinsX();
 394   for (Int_t ibin=1; ibin <= nbins; ibin++) AliInfo (Form("ibin = %d nEvProc = %g",
 395                                                           ibin,fHNProc->GetBinContent(ibin)));
 396   TFile *fout = new TFile("AliGenMUONLMR_histos.root","recreate");
 397   fHMultMu->Write();
 398   fHNProc->Write();
 399   fout->Close();
 400 }
 401
 402 //-----------------------------------------------------------
 403
 404 Double_t AliGenMUONLMR::YDistr(Double_t *px, Double_t *par){
 405   // function for rapidity distribution: plateau at par[0] +
 406   // gaussian tails centered at par[1] and with par[2]=sigma
 407   Double_t ylab = px[0];
 408   Double_t y0 = par[3]; // center of mass rapidity
 409   Double_t func = 0;
 410   if (ylab<y0+par[1] && ylab>y0-par[1]) func = par[0];
 411   else if (ylab>y0+par[1]) {
 412     Double_t z = (ylab-(par[1]+y0) )/(par[2]);
 413     func = par[0] * TMath::Exp(-0.5 * z * z);
 414   }
 415   else {
 416     Double_t z = (ylab-(-par[1]+y0) )/(par[2]);
 417     func = par[0] * TMath::Exp(-0.5 * z * z);
 418   }
 419   return func;
 420 }
 421
 422 //-----------------------------------------------------------
 423
 424 Double_t AliGenMUONLMR::PtDistr(Double_t *px, Double_t *par){
 425   // pt distribution: power law
 426   Double_t x = px[0];
 427   Double_t func = par[0] * x / TMath::Power((1+(x/par[1])*(x/par[1])),par[2]);
 428   return func;
 429 }
 430
 431 //-----------------------------------------------------------
 432
 433 void AliGenMUONLMR::Generate() {
 434   //
 435   // generate the low mass resonances and their decays according to
 436   // the multiplicity parameterized by pythia and BR from PDG
 437   // rapidity distributions parametrized from pythia
 438   // pt distributions from data (or pythia for etaprime)
 439   //
 440   Double_t pxPushed[100], pyPushed[100], pzPushed[100], ePushed[100];
 441   Int_t nmuons = -1, npartPushed = 0, pdgPushed[100];
 442   Double_t polar[3]= {0,0,0};  // Polarisation of the parent particle (for GEANT tracking)
 443   Double_t origin0[3];         // Origin of the generated parent particle (for GEANT tracking)
 444   // Calculating vertex position per event
 445   for (Int_t j=0;j<3;j++) origin0[j]=fOrigin[j];
 446   if(fVertexSmear==kPerEvent) {
 447     Vertex();
 448     for (Int_t j=0;j<3;j++) origin0[j]=fVertex[j];
 449   }
 450
 451   TParticle *mother;
 452   TDatabasePDG* pdg = TDatabasePDG::Instance();
 453
 454   Double_t pt, y, phi, mass, px, py, pz, ene, mt;
 455
 456   const Int_t nproc = 9;
 457   Int_t idRes[nproc] = {kEtaLMR, kEtaLMR, kRhoLMR, kOmegaLMR, kOmegaLMR, kPhiLMR, kEtaPrimeLMR, kPionLMR, kKaonLMR};
 458   Double_t BR[nproc] = {5.8e-6, 3.1e-4, 4.55e-5, 7.28e-5, 1.3e-4, 2.86e-4, 1.04e-4, 1, 0.6344};
 459   //  Double_t BR[nproc] = {1, 1, 1, 1, 1, 1, 1, 1, 1};
 460   Int_t idDec[nproc] = {0, 1, 0, 0, 1, 0, 1, 2, 2}; // 0:2body, 1:Dalitz, 2:pi/K
 461   Int_t mult[nproc] = {0,0,0,0,0,0,0,0,0};
 462
 463   while (nmuons < fNMuMin) {
 464
 465     nmuons = 0;
 466     npartPushed = 0;
 467     for (Int_t iproc=0; iproc<nproc; iproc++) {
 468       if (fGenSingleProc == -1) {
 469         mult[iproc] = Int_t(fMult[idRes[iproc]]->GetRandom()*fScaleMult[idRes[iproc]]);
 470       }
 471       else {
 472         if (iproc==fGenSingleProc) {
 473           mult[iproc] = 1;
 474           BR[iproc] = 1;
 475         }
 476         else {
 477           mult[iproc] = 0;
 478           BR[iproc] = 0;
 479         }
 480       }
 481     }
 482
 483     if (fGenSingleProc == -1) {
 484       mult[1] = mult[0];
 485       mult[4] = mult[3];
 486     }
 487
 488     for (Int_t iproc = 0; iproc < nproc; iproc++) {
 489       //       printf ("Multiplicity for process %d is %d\n",iproc,mult[iproc]);
 490       for (Int_t imult=0; imult<mult[iproc]; imult++) {
 491         if (gRandom->Rndm() < BR[iproc]) {
 492           fHNProc->Fill(iproc);
 493           Int_t ipart = idRes[iproc];
 494           pt  = fPt[ipart]->GetRandom();
 495           y   = fY[ipart]->GetRandom();
 496           phi = gRandom->Rndm() * 2 * TMath::Pi();
 497           mass = pdg->GetParticle(fPDG[ipart])->Mass();
 498           px  = pt * TMath::Cos(phi);
 499           py  = pt * TMath::Sin(phi);
 500           mt  = TMath::Sqrt(pt * pt + mass * mass);
 501           pz  = mt * TMath::SinH(y);
 502           ene = mt * TMath::CosH(y);
 503
 504           mother = fParticle[ipart];
 505           mother->SetMomentum(px,py,pz,ene);
 506           mother->SetCalcMass(mass);
 507           if (!KinematicSelection(mother,0)) continue;
 508
 509           Bool_t hasDecayed = kTRUE;
 510           if (idDec[iproc] == 0) Decay2Body(mother);
 511           else if (idDec[iproc] == 1) DalitzDecay(mother);
 512           else DecayPiK(mother,hasDecayed);
 513           if (!hasDecayed) continue;
 514           Bool_t isMu0Acc = KinematicSelection(fMu[0],1);
 515           Bool_t isMu1Acc = KinematicSelection(fMu[1],1);
 516           Bool_t isMuFromPiKAcc = kTRUE;
 517
 518           if (idDec[iproc] == 2) isMuFromPiKAcc = (mother->GetPdgCode()>0) ? isMu0Acc : isMu1Acc;
 519           // mother
 520           if ((idDec[iproc]  < 2 && (isMu0Acc || isMu1Acc)) ||
 521               (idDec[iproc] == 2 && isMuFromPiKAcc)) {
 522             pdgPushed[npartPushed] = mother->GetPdgCode();
 523             pxPushed[npartPushed] = mother->Px();
 524             pyPushed[npartPushed] = mother->Py();
 525             pzPushed[npartPushed] = mother->Pz();
 526             ePushed[npartPushed] = mother->Energy();
 527             npartPushed++;
 528             if (isMu0Acc && (idDec[iproc] < 2 || mother->GetPdgCode() > 0)) {
 529               pdgPushed[npartPushed] = fMu[0]->GetPdgCode();
 530               pxPushed[npartPushed] = fMu[0]->Px();
 531               pyPushed[npartPushed] = fMu[0]->Py();
 532               pzPushed[npartPushed] = fMu[0]->Pz();
 533               ePushed[npartPushed] = fMu[0]->Energy();
 534               npartPushed++;
 535               nmuons++;
 536             }
 537
 538             if (isMu1Acc && (idDec[iproc] < 2 || mother->GetPdgCode() < 0)) {
 539               pdgPushed[npartPushed] = fMu[1]->GetPdgCode();
 540               pxPushed[npartPushed] = fMu[1]->Px();
 541               pyPushed[npartPushed] = fMu[1]->Py();
 542               pzPushed[npartPushed] = fMu[1]->Pz();
 543               ePushed[npartPushed] = fMu[1]->Energy();
 544               npartPushed++;
 545               nmuons++;
 546             }
 547           }
 548         } // end if BR
 549       } // end loop on multiplicity
 550     }  // end loop on process
 551     fHMultMu->Fill(nmuons);
 552   } // keep on generating until at least a muon is created in the event
 553
 554   Int_t ntmother = 0, ntchild =0;
 555   for (Int_t ipart = 0; ipart < npartPushed; ipart++) {
 556     if (TMath::Abs(pdgPushed[ipart]) != 13) { // particle is not a muon, hence it's a mother
 557       PushTrack(0,-1,pdgPushed[ipart],
 558                 pxPushed[ipart],pyPushed[ipart],pzPushed[ipart],ePushed[ipart],
 559                 origin0[0],origin0[1],origin0[2],0.,
 560                 polar[0],polar[1],polar[2],
 561                 kPPrimary,ntmother,1,11);
 562       KeepTrack(ntmother);
 563     }
 564     else {
 565       PushTrack(1,ntmother,pdgPushed[ipart],
 566                 pxPushed[ipart],pyPushed[ipart],pzPushed[ipart],ePushed[ipart],
 567                 origin0[0],origin0[1],origin0[2],0.,
 568                 polar[0],polar[1],polar[2],
 569                 kPDecay,ntchild,1,1);
 570       KeepTrack(ntchild);
 571     }
 572   }
 573   SetHighWaterMark(ntchild);
 574   AliGenEventHeader* header = new AliGenEventHeader("LMR");
 575   header->SetPrimaryVertex(fVertex);
 576   header->SetNProduced(fNprimaries);
 577   AddHeader(header);
 578 }
 579
 580 //------------------------------------------------------------------
 581
 582 void AliGenMUONLMR::Decay2Body(TParticle *mother){
 583   // performs decay in two muons of the low mass resonances
 584   Double_t md1 = fMu[0]->GetMass();
 585   Int_t pdg = mother->GetPdgCode();
 586   Double_t mres =0;
 587   // if mother is a rho, extract the mass from its line shape
 588   // otherwise consider the resonance mass
 589   if (pdg == 113) mres = fRhoLineShape->GetRandom();
 590   else mres = mother->GetCalcMass();
 591   //  while (mres < md1 + md2) mres =  fDsigmaDm[res]->GetRandom();
 592   // energies and momenta in rest frame
 593   Double_t e1 = mres / 2.;
 594   Double_t p1 = TMath::Sqrt((e1 + md1)*(e1 - md1));
 595   // orientation in decaying particle rest frame
 596   Double_t costheta = gRandom->Rndm() * 2 - 1;
 597   Double_t sintheta = TMath::Sqrt((1. + costheta)*(1. - costheta));
 598   Double_t phi      = 2. * TMath::Pi() * gRandom->Rndm();
 599   Double_t px1      = p1 * sintheta * TMath::Cos(phi);
 600   Double_t py1      = p1 * sintheta * TMath::Sin(phi);
 601   Double_t pz1      = p1 * costheta;
 602
 603   // boost muons into lab frame
 604
 605   TLorentzVector vmother, v1, v2;
 606   //  TLorentzVector boosted1, boosted2;
 607   vmother.SetPxPyPzE(mother->Px(),mother->Py(),mother->Pz(),mother->Energy());
 608   v1.SetPxPyPzE(px1,py1,pz1,e1);
 609   v2.SetPxPyPzE(-px1,-py1,-pz1,e1);
 610
 611   TVector3 betaParent = (1./vmother.E())*vmother.Vect(); // beta = p/E
 612   v1.Boost(betaParent);
 613   v2.Boost(betaParent);
 614
 615   //   TLorentzVector vtot = v1 + v2;
 616   //   printf ("mother: %g   %g    %g     %g\n",vmother.Px(), vmother.Py(), vmother.Pz(), vmother.E());
 617   //   printf ("vtot  : %g   %g    %g     %g\n",vtot.Px(), vtot.Py(), vtot.Pz(), vtot.E());
 618
 619   fMu[0]->SetMomentum(v1.Px(),v1.Py(),v1.Pz(),v1.E());
 620   fMu[1]->SetMomentum(v2.Px(),v2.Py(),v2.Pz(),v2.E());
 621 }
 622
 623 //------------------------------------------------------------------
 624
 625 void AliGenMUONLMR::DecayPiK(TParticle *mother, Bool_t &hasDecayed){
 626   // performs decays of pions and kaons
 627   Double_t md1 = fMu[0]->GetMass();
 628   // extract the mass from the resonance's line shape
 629   Double_t mres = mother->GetMass();
 630   // choose the pi/k sign, assuming 50% probabilities for both signs
 631   Int_t sign = (gRandom->Rndm() > 0.5) ? 1 : -1;
 632   mother->SetPdgCode(sign * TMath::Abs(mother->GetPdgCode()));
 633
 634   // energies and momenta in rest frame
 635   Double_t e1 = (mres*mres + md1*md1)/(2*mres);
 636   Double_t p1 = TMath::Sqrt((e1 + md1)*(e1 - md1));
 637   // orientation in decaying particle rest frame
 638   Double_t costheta = gRandom->Rndm() * 2 - 1;
 639   Double_t sintheta = TMath::Sqrt((1. + costheta)*(1. - costheta));
 640   Double_t phi      = 2. * TMath::Pi() * gRandom->Rndm();
 641   Double_t px1      = p1 * sintheta * TMath::Cos(phi);
 642   Double_t py1      = p1 * sintheta * TMath::Sin(phi);
 643   Double_t pz1      = p1 * costheta;
 644
 645   // boost muons into lab frame
 646   TLorentzVector vmother, v1;
 647   vmother.SetPxPyPzE(mother->Px(),mother->Py(),mother->Pz(),mother->Energy());
 648   v1.SetPxPyPzE(px1,py1,pz1,e1);
 649
 650   TVector3 betaParent = (1./vmother.E())*vmother.Vect(); // beta = p/E
 651   v1.Boost(betaParent);
 652   if (mother->GetPdgCode()>0) fMu[0]->SetMomentum(v1.Px(),v1.Py(),v1.Pz(),v1.E());
 653   else fMu[1]->SetMomentum(v1.Px(),v1.Py(),v1.Pz(),v1.E());
 654
 655   Int_t idmother = 0;
 656   if (TMath::Abs(mother->GetPdgCode())== 211) idmother = 0;
 657   if (TMath::Abs(mother->GetPdgCode())== 321) idmother = 1;
 658   Double_t gammaRes = mother->Energy()/mres;
 659   Double_t zResCM = fDecay[idmother]->GetRandom();
 660   Double_t zResLab = gammaRes*zResCM;
 661   if(zResLab > 0.938) hasDecayed = 0; // 0.938: distance from IP to absorber + lambda_i
 662   else hasDecayed = 1;
 663 }
 664
 665 //-------------------------------------------------------------------
 666
 667 void AliGenMUONLMR::DalitzDecay(TParticle *mother){
 668   //
 669   // perform dalitz decays of eta, omega and etaprime
 670   //
 671   //in the rest frame of the virtual photon:
 672   Double_t mres = mother->GetCalcMass();
 673   Double_t mumass  = fMu[0]->GetMass();
 674   Double_t md3  = 0;  // unless differently specified, third particle is a photon
 675   if (mother->GetPdgCode() == 223) md3 = 0.134977; // if mother is an omega, third particle is a pi0
 676   Int_t index = 0;
 677   if (mother->GetPdgCode() == 221) index = 0;  // eta
 678   else if (mother->GetPdgCode() == 223) index = 1; // omega
 679   else if (mother->GetPdgCode() == 331) index = 2; // etaPrime
 680   Int_t flag = 0;
 681   Double_t xd=0, mvirt2=0;
 682   Double_t countIt = 0;
 683   while (flag==0) {
 684     xd       = fDalitz[index]->GetRandom();
 685     mvirt2   = xd * mres * mres;   // mass of virtual photon
 686     // check kinematics
 687     if (mres - md3 > TMath::Sqrt(mvirt2) && TMath::Sqrt(mvirt2)/2. > mumass) flag=1;
 688     if (++countIt>1E11) {
 689       mvirt2 =  mres * mres * 0.998;
 690       break;
 691     }
 692   }
 693
 694   //
 695   //        Generate muons in virtual photon rest frame.
 696   //        z axis is the virt. photon direction (before boost)
 697   //
 698
 699   Double_t e1 = TMath::Sqrt(mvirt2)/2.; // energy of mu1 in the virtual photon frame
 700   Double_t psquare = (e1 + mumass)*(e1 - mumass);
 701   if (psquare<0) {
 702     AliError(Form("sqrt of psquare = %f put to 0\n",psquare));
 703     psquare = 0;
 704   }
 705   Double_t p1 = TMath::Sqrt(psquare);
 706   //theta angle between the pos. muon and the virtual photon
 707   Double_t costheta = fCosTheta->GetRandom();
 708   if (costheta>1)  costheta = 1;
 709   if (costheta<-1) costheta = -1;
 710   Double_t sintheta = TMath::Sqrt((1. + costheta)*(1. - costheta));
 711   Double_t phi      = 2 * TMath::Pi() * gRandom->Rndm();
 712   Double_t sinphi   = TMath::Sin(phi);
 713   Double_t cosphi   = TMath::Cos(phi);
 714
 715   // fill 4-vectors of leptons in the virtual photon frame
 716
 717   Double_t px1 = p1*sintheta*cosphi;
 718   Double_t py1 = p1*sintheta*sinphi;
 719   Double_t pz1 = p1*costheta;
 720   Double_t px2 = -p1*sintheta*cosphi;
 721   Double_t py2 = -p1*sintheta*sinphi;
 722   Double_t pz2 = -p1*costheta;
 723   Double_t e2  = e1;
 724
 725   fMu[0]->SetMomentum(px1,py1,pz1,e1);
 726   fMu[1]->SetMomentum(px2,py2,pz2,e2);
 727
 728   // calculate components of non-dilepton in CMS of parent resonance
 729
 730   Double_t e3 = (mres * mres + md3 * md3 - mvirt2) / (2.*mres);
 731   Double_t psquare3 = (e3 + md3)*(e3 - md3);
 732   if (psquare3<0) {
 733     AliError(Form("Sqrt of psquare3 = %f put to 0\n",psquare3));
 734     psquare3 = 0;
 735   }
 736   Double_t p3 = TMath::Sqrt(psquare3);
 737   Double_t costheta2 = 2.* gRandom->Rndm() - 1.;   // angle between virtual photon and resonance
 738   if (costheta2>1)  costheta2 = 1;
 739   if (costheta2<-1) costheta2 = -1;
 740   Double_t sintheta2 = TMath::Sqrt((1. + costheta2)*(1. - costheta2));
 741   Double_t phi2      = 2 * TMath::Pi() * gRandom->Rndm();
 742   Double_t sinphi2   = TMath::Sin(phi2);
 743   Double_t cosphi2   = TMath::Cos(phi2);
 744   Double_t px3 = p3*sintheta2*cosphi2;
 745   Double_t py3 = p3*sintheta2*sinphi2;
 746   Double_t pz3 = p3*costheta2;
 747   TLorentzVector v3(px3,py3,pz3,e3);
 748
 749   sintheta2 = -sintheta2;
 750   cosphi2   = -cosphi2;
 751   sinphi2   = -sinphi2;
 752
 753   Double_t px1new = px1*costheta2*cosphi2 - py1*sinphi2 + pz1*sintheta2*cosphi2;
 754   Double_t py1new = px1*costheta2*sinphi2 + py1*cosphi2 + pz1*sintheta2*sinphi2;
 755   Double_t pz1new =-px1*sintheta2                       + pz1*costheta2;
 756   Double_t px2new = px2*costheta2*cosphi2 - py2*sinphi2 + pz2*sintheta2*cosphi2;
 757   Double_t py2new = px2*costheta2*sinphi2 + py2*cosphi2 + pz2*sintheta2*sinphi2;
 758   Double_t pz2new =-px2*sintheta2                       + pz2*costheta2;
 759
 760   fMu[0]->SetMomentum(px1new,py1new,pz1new,e1);
 761   fMu[1]->SetMomentum(px2new,py2new,pz2new,e2);
 762
 763   Double_t evirt = mres - e3;
 764   Double_t pxvirt = -px3;
 765   Double_t pyvirt = -py3;
 766   Double_t pzvirt = -pz3;
 767   TLorentzVector vvirt(pxvirt,pyvirt,pzvirt,evirt);
 768   TVector3 betaVirt = (1./evirt) * vvirt.Vect(); // virtual photon beta in res frame
 769
 770   TLorentzVector v1(px1,py1,pz1,e1);
 771   TLorentzVector v2(px2,py2,pz2,e2);
 772
 773   // boost the muons in the frame where the resonance is at rest
 774
 775   v1.Boost(betaVirt);
 776   v2.Boost(betaVirt);
 777
 778   // boost muons and third particle in lab frame
 779
 780   TLorentzVector vmother(mother->Px(), mother->Py(), mother->Pz(), mother->Energy());
 781   TVector3 resBetaLab = (1./vmother.E())*vmother.Vect(); // eta beta in lab frame
 782   v1.Boost(resBetaLab);
 783   v2.Boost(resBetaLab);
 784   v3.Boost(resBetaLab);
 785   vvirt.Boost(resBetaLab);
 786
 787   fMu[0]->SetMomentum(v1.Px(),v1.Py(),v1.Pz(),v1.E());
 788   fMu[1]->SetMomentum(v2.Px(),v2.Py(),v2.Pz(),v2.E());
 789   //   part3->SetMomentum(v3.Px(),v3.Py(),v3.Pz(),v3.E());
 790
 791   //   TLorentzVector vtot = v1 + v2 + v3;
 792   //   TLorentzVector vdimu = v1 + v2;
 793   //   printf ("mother: %g   %g    %g     %g\n",vmother.Px(), vmother.Py(), vmother.Pz(), vmother.E());
 794   //   printf ("vtot  : %g   %g    %g     %g\n",vtot.Px(), vtot.Py(), vtot.Pz(), vtot.E());
 795   //   printf ("vvirt : %g   %g    %g     %g\n",vvirt.Px(), vvirt.Py(), vvirt.Pz(), vvirt.E());
 796   //   printf ("vdimu : %g   %g    %g     %g\n",vdimu.Px(), vdimu.Py(), vdimu.Pz(), vdimu.E());
 797
 798 }
 799
 800 //------------------------------------------------------------------
 801
 802 Double_t AliGenMUONLMR::FormFactor(Double_t q2, Int_t decay){
 803   //  Calculates the form factor for Dalitz decays A->B+l+l
 804   //  Returns: |F(q^2)|^2
 805   //
 806   //  References: L.G. Landsberg, Physics Reports 128 No.6 (1985) 301-376.
 807
 808   Double_t ff2, mass2;
 809   Double_t n2, n4, m2;
 810   // Lepton-G
 811
 812   Double_t lambda2inv = 0;
 813   switch (decay) {
 814   case 0:   // eta -> mu mu gamma
 815     // eta   -> l+ l- gamma: pole approximation
 816     lambda2inv = 1.95;
 817     mass2 = fParticle[kEtaLMR]->GetMass() * fParticle[kEtaLMR]->GetMass();
 818     if (q2 < mass2) ff2 = TMath::Power(1./(1.-lambda2inv*q2),2);
 819     else ff2 = 0;
 820     break;
 821   case 1:   // omega -> mu mu pi0
 822     // omega -> l+ l- pi0: pole approximation
 823     mass2 = fParticle[kOmegaLMR]->GetMass() * fParticle[kOmegaLMR]->GetMass();
 824     lambda2inv = 2.26;
 825     if (q2 < mass2) ff2 = TMath::Power(1./(1.-lambda2inv*q2),2);
 826     else ff2 = 0;
 827     break;
 828   case 2:   // etaPrime -> mu mu gamma
 829     mass2 = fParticle[kEtaPrimeLMR]->GetMass() * fParticle[kEtaPrimeLMR]->GetMass();
 830     // eta'  -> l+ l- gamma: Breit-Wigner fitted to data
 831     n2 = 0.764 * 0.764;
 832     n4 = n2 * n2;
 833     m2 = 0.1020 * 0.1020;
 834     if (q2 < mass2) ff2 = n4 / (TMath::Power(n2-q2,2) + m2 * n2);
 835     else ff2 = 0;
 836     break;
 837   default:
 838     AliError ("FormFactor: Decay not found");
 839     return 0;
 840     break;
 841   }
 842   return ff2;
 843 }
 844
 845 //____________________________________________________________
 846
 847 Double_t AliGenMUONLMR::RhoLineShapeNew(Double_t *x, Double_t */*para*/){
 848   //new parameterization implemented by Hiroyuki Sako (GSI)
 849   Double_t mass = *x;
 850   double r, GammaTot;
 851   Double_t mRho    = TDatabasePDG::Instance()->GetParticle("rho0")->Mass();
 852   Double_t mPi     = TDatabasePDG::Instance()->GetParticle("pi0")->Mass();
 853   Double_t mMu     = TDatabasePDG::Instance()->GetParticle("mu-")->Mass();
 854   Double_t Gamma0  = TDatabasePDG::Instance()->GetParticle("rho0")->Width();
 855
 856   const double Norm = 0.0744416*1.01;
 857
 858   // 0.0744416 at m = 0.72297
 859   // is the max number with Norm=1 (for rho)
 860
 861   double mThreshold = 2.*mPi;
 862
 863   const double T = 0.170; // Assumption of pi+ temperature [GeV/c^2]
 864   //const double T = 0.11; // Taken from fit to pi+ temperature [GeV/c^2]
 865   // with Reference: LEBC-EHS collab., Z. Phys. C 50 (1991) 405
 866
 867   if (mass < mThreshold) {
 868     r = 0.;
 869     return r;
 870   }
 871
 872   double k = sqrt(0.25*mass*mass-(mThreshold/2)*(mThreshold/2));
 873   double k0 = sqrt(0.25*mRho*mRho-(mThreshold/2)*(mThreshold/2));
 874
 875   GammaTot = (k/k0)*(k/k0)*(k/k0)*(mRho/mass)*(mRho/mass)*Gamma0;
 876
 877   double FormFactor2 = 1/((mass*mass-mRho*mRho)*(mass*mass-mRho*mRho)+
 878                           mass*mass*GammaTot*GammaTot);
 879
 880   r = pow(mass,1.5)*pow((1-mThreshold*mThreshold/(mass*mass)),1.5)*
 881     ((mass*mass+2*mMu*mMu)/(mass*mass))*(pow((mass*mass-4*mMu*mMu),0.5)/mass)*FormFactor2
 882     *exp(-mass/T)/Norm;
 883
 884   return r;
 885 }
 886