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