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