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