EVGEN/AliGenMUONLMR.cxx

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