X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EVGEN%2FAliGenMUONLMR.cxx;h=5408ec8f500fbeb66be11a3787428b88ffebb9c3;hb=dd76cfd5dcacd53211f68c82c6d3c406422a518a;hp=602991a8ee4bff966f2c0a1d31d095a8cf1cd848;hpb=4b4ed190878fb59b27b03cfb0b8efe356d77352e;p=u%2Fmrichter%2FAliRoot.git diff --git a/EVGEN/AliGenMUONLMR.cxx b/EVGEN/AliGenMUONLMR.cxx index 602991a8ee4..5408ec8f500 100644 --- a/EVGEN/AliGenMUONLMR.cxx +++ b/EVGEN/AliGenMUONLMR.cxx @@ -3,57 +3,215 @@ #include #include #include "AliConst.h" -#include "AliLog.h" #include "AliGenMUONLMR.h" #include "AliMC.h" #include "AliRun.h" +#include "AliLog.h" #include "AliGenEventHeader.h" ClassImp(AliGenMUONLMR) + AliGenMUONLMR::AliGenMUONLMR () : AliGenMC(), + fNMuMin(2), + fCMSEnergy(kNCMSEnergies), + fGenSingleProc(-1), + fYCM(0), + fCosTheta(0x0), + fRhoLineShape(0x0), + fHMultMu(0x0), + fHNProc(0x0) { + // + // default constructor + // + for (int i=0; iALICE_2.76 = ALICE_7 * PYTHIA_2.76 / PYTHIA_7 + AliInfo ("Using pp parameterization at 2.76 TeV\n"); + Double_t yparam2760[7][9] = {{1,0.8251,3.657,0,0,0,0,0,0}, // pions from pythia + {1,1.83,2.698,0,0,0,0,0,0}, // kaons from pythia + {1,0.011,3.474,0,0,0,0,0,0}, // eta from pythia + {1,-0.01,3.409,0,0,0,0,0,0}, // rho from pythia + {1,-0.037,3.294,0,0,0,0,0,0}, // omega from pythia + {1,-0.016,2.717,0,0,0,0,0,0}, // phi from pythia + {1,-0.010,3.312,0,0,0,0,0,0}}; // eta prime from pythia + + Double_t ptparam2760[7][9] = {{1,0.1665,8.878,0,0,0,0,0,0}, // pions from Pythia + {1,0.1657,8.591,0,0,0,0,0,0}, // kaons from Pythia + {1,0.641,2.62,0,0,0,0,0,0}, // eta from ALICE 7 TeV + {1,1.3551,3.16,0,0,0,0,0,0}, // rho with scaled + {1,1.3551,3.16,0,0,0,0,0,0}, // omega with scaled + {1,1.0811,2.74,0,0,0,0,0,0}, // phi with scaled + {1,0.72,2.5,0,0,0,0,0,0}}; // etaPrime from ALICE 7 TeV + + Double_t nparam2760[7][9] = {{9752,-2.693,3.023,9.5e9,-84.68,16.75,-14.06,635.3,-423.2}, // pions + {1.e5, 1.538, 0,0,0,0,0,0,0}, // kaons + {1.e4, 0.351, 0,0,0,0,0,0,0}, // eta + {1.e4, 0.2471, 0,0,0,0,0,0,0}, // rho + {1.e4, 0.2583, 0,0,0,0,0,0,0}, // omega + {1.e5, 1.393, 0,0,0,0,0,0,0}, // phi + {1.e4, 0.9005, 0,0,0,0,0,0,0}}; // etaPrime + + for (Int_t i=0; iSetParameter(j,nparam[i][j]); } - + fPt[i] = new TF1(fptname[i],AliGenMUONLMR::PtDistr,0,20,3); fPt[i]->SetParameters(ptparam[i][0], ptparam[i][1], ptparam[i][2]); - fY[i] = new TF1(fyname[i],AliGenMUONLMR::YDistr,-10,10,3); - fY[i]->SetParameters(yparam[i][0], yparam[i][1], yparam[i][2]); + fY[i] = new TF1(fyname[i],AliGenMUONLMR::YDistr,-10,10,4); + fY[i]->SetParameters(yparam[i][0], yparam[i][1], yparam[i][2],fYCM); } - + for(Int_t i = 0; i<2; i++){ fDecay[i] = new TF1(fdname[i],"exp(-x/[0])",0,150); fDecay[i]->SetParameter(0,ctau[i]); } - + for (Int_t ipart = 0; ipart < fgkNpart; ipart++) { fParticle[ipart] = new TParticle(); fParticle[ipart]->SetPdgCode(fPDG[ipart]); } - + TDatabasePDG *pdgdb = TDatabasePDG::Instance(); Double_t mumass = pdgdb->GetParticle(13)->Mass(); fMu[0] = new TParticle(); @@ -89,24 +247,24 @@ AliGenMUONLMR::AliGenMUONLMR () : AliGenMC(), fNMuMin(2), fGenSingleProc(-1), fC fMu[1] = new TParticle(); fMu[1]->SetPdgCode(13); fMu[1]->SetCalcMass(mumass); - + // function for polarized theta distributions fCosTheta = new TF1 ("fCosTheta","1+[0]*x*x",-1,1); fCosTheta->SetParameter(0,1); - + // Dalitz decays Int_t nbins = 1000; Double_t xmin = 0, xmax = 2; fDalitz[0] = new TH1F("hDalitzEta","",nbins,xmin,xmax); fDalitz[1] = new TH1F("hDalitzOmega","",nbins,xmin,xmax); fDalitz[2] = new TH1F("hDalitzEtaPrime","",nbins,xmin,xmax); - + Double_t meta = pdgdb->GetParticle("eta")->Mass(); Double_t momega = pdgdb->GetParticle("omega")->Mass(); Double_t metaPrime = pdgdb->GetParticle("eta'")->Mass(); Double_t mpi0 = pdgdb->GetParticle("pi0")->Mass(); Double_t md3 = 0, mres = 0; - + for (Int_t index = 0; index < 3; index++) { if (index == 0) { mres = meta; @@ -122,11 +280,11 @@ AliGenMUONLMR::AliGenMUONLMR () : AliGenMC(), fNMuMin(2), fGenSingleProc(-1), fC } Double_t delta = md3 * md3 / (mres * mres); Double_t epsilon = mumass * mumass / (mres * mres); - Int_t nbins0 = fDalitz[index]->GetNbinsX(); - Double_t xmin0 = fDalitz[index]->GetXaxis()->GetXmin(); + nbins = fDalitz[index]->GetNbinsX(); + xmin = fDalitz[index]->GetXaxis()->GetXmin(); Double_t deltax = fDalitz[index]->GetBinWidth(1); - Double_t xd = xmin0 - deltax/2.; - for (Int_t ibin = 0; ibin< nbins0; ibin++) { + Double_t xd = xmin - deltax/2.; + for (Int_t ibin = 0; ibin< nbins; ibin++) { Double_t dalval = 0; xd += deltax; if (xd > 4. *epsilon) { @@ -141,12 +299,68 @@ AliGenMUONLMR::AliGenMUONLMR () : AliGenMC(), fNMuMin(2), fGenSingleProc(-1), fC } } } - + fRhoLineShape = new TF1("fRhoLineShape",RhoLineShapeNew,0,2,2); fHMultMu = new TH1D("fHMultMu","Muon multiplicity",20,-0.5,19.5); fHNProc = new TH1D("fHNProc","Number of gen. evts. per process in 4 pi",9,-0.5,8.5); + +} +//----------------------------------------------------------- + +AliGenMUONLMR::AliGenMUONLMR (AliGenMUONLMR &gen) : AliGenMC(), + fNMuMin(gen.fNMuMin), + fCMSEnergy(gen.fCMSEnergy), + fGenSingleProc(gen.fGenSingleProc), + fYCM(gen.fYCM), + fCosTheta(gen.fCosTheta), + fRhoLineShape(gen.fRhoLineShape), + fHMultMu(gen.fHMultMu), + fHNProc(gen.fHNProc) { + for (Int_t i=0; i < fgkNpart; i++) { + fPDG[i] = gen.fPDG[i]; + fScaleMult[i] = gen.fScaleMult[i]; + fPt[i] = (TF1*) gen.fPt[i]->Clone(); + fY[i] = (TF1*) gen.fY[i]->Clone(); + fMult[i] = (TF1*) gen.fMult[i]->Clone(); + fParticle[i] = (TParticle*) gen.fParticle[i]->Clone(); + } + + for(Int_t i = 0; i<2; i++) fDecay[i] = (TF1*) gen.fDecay[i]->Clone(); + for(Int_t i = 0; i<3; i++) fDalitz[i] = (TH1F*) gen.fDalitz[i]->Clone(); + for(Int_t i = 0; i<2; i++) fMu[i] = (TParticle*) gen.fMu[i]->Clone(); } +//----------------------------------------------------------- + +AliGenMUONLMR& AliGenMUONLMR::operator=(const AliGenMUONLMR &gen) { + // Assignment operator + if (this!=&gen) { + fNMuMin = gen.fNMuMin; + fCMSEnergy = gen.fCMSEnergy; + fGenSingleProc = gen.fGenSingleProc; + fYCM = gen.fYCM; + fCosTheta = (TF1*) gen.fCosTheta->Clone(); + fRhoLineShape = (TF1*) gen.fRhoLineShape->Clone(); + fHMultMu = (TH1D*) gen.fHMultMu->Clone(); + fHNProc = (TH1D*) gen.fHNProc->Clone(); + + for (Int_t i=0; i < fgkNpart; i++) { + fPDG[i] = gen.fPDG[i]; + fScaleMult[i] = gen.fScaleMult[i]; + fPt[i] = (TF1*) gen.fPt[i]->Clone(); + fY[i] = (TF1*) gen.fY[i]->Clone(); + fMult[i] = (TF1*) gen.fMult[i]->Clone(); + fParticle[i] = (TParticle*) gen.fParticle[i]->Clone(); + } + + for(Int_t i = 0; i<2; i++) fDecay[i] = (TF1*) gen.fDecay[i]->Clone(); + for(Int_t i = 0; i<3; i++) fDalitz[i] = (TH1F*) gen.fDalitz[i]->Clone(); + for(Int_t i = 0; i<2; i++) fMu[i] = (TParticle*) gen.fMu[i]->Clone(); + } + return *this; +} + + //----------------------------------------------------------- AliGenMUONLMR::~AliGenMUONLMR() @@ -177,9 +391,9 @@ AliGenMUONLMR::~AliGenMUONLMR() void AliGenMUONLMR::FinishRun(){ // save some histograms to an output file Int_t nbins = fHNProc->GetNbinsX(); - for (Int_t ibin=1; ibin <= nbins; ibin++) printf ("ibin = %d nEvProc = %g\n", - ibin,fHNProc->GetBinContent(ibin)); - TFile *fout = new TFile("AliGenMUONLMR_histos.root","recreate"); + for (Int_t ibin=1; ibin <= nbins; ibin++) AliInfo (Form("ibin = %d nEvProc = %g", + ibin,fHNProc->GetBinContent(ibin))); + TFile *fout = new TFile("AliGenMUONLMR_histos.root","recreate"); fHMultMu->Write(); fHNProc->Write(); fout->Close(); @@ -187,14 +401,19 @@ void AliGenMUONLMR::FinishRun(){ //----------------------------------------------------------- -Double_t AliGenMUONLMR::YDistr(const Double_t *px, const Double_t *par){ +Double_t AliGenMUONLMR::YDistr(Double_t *px, Double_t *par){ // function for rapidity distribution: plateau at par[0] + // gaussian tails centered at par[1] and with par[2]=sigma - Double_t x = TMath::Abs(px[0]); + Double_t ylab = px[0]; + Double_t y0 = par[3]; // center of mass rapidity Double_t func = 0; - if (xy0-par[1]) func = par[0]; + else if (ylab>y0+par[1]) { + Double_t z = (ylab-(par[1]+y0) )/(par[2]); + func = par[0] * TMath::Exp(-0.5 * z * z); + } + else { + Double_t z = (ylab-(-par[1]+y0) )/(par[2]); func = par[0] * TMath::Exp(-0.5 * z * z); } return func; @@ -202,7 +421,7 @@ Double_t AliGenMUONLMR::YDistr(const Double_t *px, const Double_t *par){ //----------------------------------------------------------- -Double_t AliGenMUONLMR::PtDistr(const Double_t *px, const Double_t *par){ +Double_t AliGenMUONLMR::PtDistr(Double_t *px, Double_t *par){ // pt distribution: power law Double_t x = px[0]; Double_t func = par[0] * x / TMath::Power((1+(x/par[1])*(x/par[1])),par[2]); @@ -222,17 +441,13 @@ void AliGenMUONLMR::Generate() { Int_t nmuons = -1, npartPushed = 0, pdgPushed[100]; Double_t polar[3]= {0,0,0}; // Polarisation of the parent particle (for GEANT tracking) Double_t origin0[3]; // Origin of the generated parent particle (for GEANT tracking) - Double_t time0; // Time0 of the generated parent particle // Calculating vertex position per event for (Int_t j=0;j<3;j++) origin0[j]=fOrigin[j]; - time0 = fTimeOrigin; if(fVertexSmear==kPerEvent) { Vertex(); for (Int_t j=0;j<3;j++) origin0[j]=fVertex[j]; - time0 = fTime; } - printf ("In Generate()\n"); TParticle *mother; TDatabasePDG* pdg = TDatabasePDG::Instance(); @@ -240,7 +455,7 @@ void AliGenMUONLMR::Generate() { const Int_t nproc = 9; Int_t idRes[nproc] = {kEtaLMR, kEtaLMR, kRhoLMR, kOmegaLMR, kOmegaLMR, kPhiLMR, kEtaPrimeLMR, kPionLMR, kKaonLMR}; - 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}; + 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}; // Double_t BR[nproc] = {1, 1, 1, 1, 1, 1, 1, 1, 1}; Int_t idDec[nproc] = {0, 1, 0, 0, 1, 0, 1, 2, 2}; // 0:2body, 1:Dalitz, 2:pi/K Int_t mult[nproc] = {0,0,0,0,0,0,0,0,0}; @@ -271,7 +486,7 @@ void AliGenMUONLMR::Generate() { } for (Int_t iproc = 0; iproc < nproc; iproc++) { -// printf ("Multiplicity for process %d is %d\n",iproc,mult[iproc]); + // printf ("Multiplicity for process %d is %d\n",iproc,mult[iproc]); for (Int_t imult=0; imultRndm() < BR[iproc]) { fHNProc->Fill(iproc); @@ -341,7 +556,7 @@ void AliGenMUONLMR::Generate() { if (TMath::Abs(pdgPushed[ipart]) != 13) { // particle is not a muon, hence it's a mother PushTrack(0,-1,pdgPushed[ipart], pxPushed[ipart],pyPushed[ipart],pzPushed[ipart],ePushed[ipart], - origin0[0],origin0[1],origin0[2],time0, + origin0[0],origin0[1],origin0[2],0., polar[0],polar[1],polar[2], kPPrimary,ntmother,1,11); KeepTrack(ntmother); @@ -349,7 +564,7 @@ void AliGenMUONLMR::Generate() { else { PushTrack(1,ntmother,pdgPushed[ipart], pxPushed[ipart],pyPushed[ipart],pzPushed[ipart],ePushed[ipart], - origin0[0],origin0[1],origin0[2],time0, + origin0[0],origin0[1],origin0[2],0., polar[0],polar[1],polar[2], kPDecay,ntchild,1,1); KeepTrack(ntchild); @@ -358,14 +573,13 @@ void AliGenMUONLMR::Generate() { SetHighWaterMark(ntchild); AliGenEventHeader* header = new AliGenEventHeader("LMR"); header->SetPrimaryVertex(fVertex); - header->SetInteractionTime(fTime); header->SetNProduced(fNprimaries); AddHeader(header); } //------------------------------------------------------------------ -void AliGenMUONLMR::Decay2Body(const TParticle *mother){ +void AliGenMUONLMR::Decay2Body(TParticle *mother){ // performs decay in two muons of the low mass resonances Double_t md1 = fMu[0]->GetMass(); Int_t pdg = mother->GetPdgCode(); @@ -398,9 +612,9 @@ void AliGenMUONLMR::Decay2Body(const TParticle *mother){ v1.Boost(betaParent); v2.Boost(betaParent); -// TLorentzVector vtot = v1 + v2; -// printf ("mother: %g %g %g %g\n",vmother.Px(), vmother.Py(), vmother.Pz(), vmother.E()); -// printf ("vtot : %g %g %g %g\n",vtot.Px(), vtot.Py(), vtot.Pz(), vtot.E()); + // TLorentzVector vtot = v1 + v2; + // printf ("mother: %g %g %g %g\n",vmother.Px(), vmother.Py(), vmother.Pz(), vmother.E()); + // printf ("vtot : %g %g %g %g\n",vtot.Px(), vtot.Py(), vtot.Pz(), vtot.E()); fMu[0]->SetMomentum(v1.Px(),v1.Py(),v1.Pz(),v1.E()); fMu[1]->SetMomentum(v2.Px(),v2.Py(),v2.Pz(),v2.E()); @@ -439,25 +653,18 @@ void AliGenMUONLMR::DecayPiK(TParticle *mother, Bool_t &hasDecayed){ else fMu[1]->SetMomentum(v1.Px(),v1.Py(),v1.Pz(),v1.E()); Int_t idmother = 0; - if (TMath::Abs(mother->GetPdgCode())== 211) { - idmother = 0; - } else if (TMath::Abs(mother->GetPdgCode())== 321) { - idmother = 1; - } else { - AliWarning("Mother particle is not a pion or kaon \n"); - } - + if (TMath::Abs(mother->GetPdgCode())== 211) idmother = 0; + if (TMath::Abs(mother->GetPdgCode())== 321) idmother = 1; Double_t gammaRes = mother->Energy()/mres; Double_t zResCM = fDecay[idmother]->GetRandom(); Double_t zResLab = gammaRes*zResCM; if(zResLab > 0.938) hasDecayed = 0; // 0.938: distance from IP to absorber + lambda_i else hasDecayed = 1; - } //------------------------------------------------------------------- -void AliGenMUONLMR::DalitzDecay(const TParticle *mother){ +void AliGenMUONLMR::DalitzDecay(TParticle *mother){ // // perform dalitz decays of eta, omega and etaprime // @@ -466,21 +673,10 @@ void AliGenMUONLMR::DalitzDecay(const TParticle *mother){ Double_t mumass = fMu[0]->GetMass(); Double_t md3 = 0; // unless differently specified, third particle is a photon if (mother->GetPdgCode() == 223) md3 = 0.134977; // if mother is an omega, third particle is a pi0 - Int_t index = 0; - if (TMath::Abs(mother->GetPdgCode())== 221) { - // eta - index = 0; - } else if (TMath::Abs(mother->GetPdgCode())== 223) { - // omega - index = 1; - } else if (mother->GetPdgCode() == 331) { - // eta' - index = 2; - } else { - AliWarning("Mother particle is not a eta, eta' or omega \n"); - } - + if (mother->GetPdgCode() == 221) index = 0; // eta + else if (mother->GetPdgCode() == 223) index = 1; // omega + else if (mother->GetPdgCode() == 331) index = 2; // etaPrime Int_t flag = 0; Double_t xd=0, mvirt2=0; Double_t countIt = 0; @@ -503,7 +699,7 @@ void AliGenMUONLMR::DalitzDecay(const TParticle *mother){ Double_t e1 = TMath::Sqrt(mvirt2)/2.; // energy of mu1 in the virtual photon frame Double_t psquare = (e1 + mumass)*(e1 - mumass); if (psquare<0) { - printf("Error in AliGenMUONLMR::DalitzDecay: sqrt of psquare = %f put to 0\n",psquare); + AliError(Form("sqrt of psquare = %f put to 0\n",psquare)); psquare = 0; } Double_t p1 = TMath::Sqrt(psquare); @@ -534,7 +730,7 @@ void AliGenMUONLMR::DalitzDecay(const TParticle *mother){ Double_t e3 = (mres * mres + md3 * md3 - mvirt2) / (2.*mres); Double_t psquare3 = (e3 + md3)*(e3 - md3); if (psquare3<0) { - printf("Error in AliGenMUONLMR::DalitzDecay: sqrt of psquare3 = %f put to 0\n",psquare3); + AliError(Form("Sqrt of psquare3 = %f put to 0\n",psquare3)); psquare3 = 0; } Double_t p3 = TMath::Sqrt(psquare3); @@ -590,14 +786,14 @@ void AliGenMUONLMR::DalitzDecay(const TParticle *mother){ fMu[0]->SetMomentum(v1.Px(),v1.Py(),v1.Pz(),v1.E()); fMu[1]->SetMomentum(v2.Px(),v2.Py(),v2.Pz(),v2.E()); -// part3->SetMomentum(v3.Px(),v3.Py(),v3.Pz(),v3.E()); + // part3->SetMomentum(v3.Px(),v3.Py(),v3.Pz(),v3.E()); -// TLorentzVector vtot = v1 + v2 + v3; -// TLorentzVector vdimu = v1 + v2; -// printf ("mother: %g %g %g %g\n",vmother.Px(), vmother.Py(), vmother.Pz(), vmother.E()); -// printf ("vtot : %g %g %g %g\n",vtot.Px(), vtot.Py(), vtot.Pz(), vtot.E()); -// printf ("vvirt : %g %g %g %g\n",vvirt.Px(), vvirt.Py(), vvirt.Pz(), vvirt.E()); -// printf ("vdimu : %g %g %g %g\n",vdimu.Px(), vdimu.Py(), vdimu.Pz(), vdimu.E()); + // TLorentzVector vtot = v1 + v2 + v3; + // TLorentzVector vdimu = v1 + v2; + // printf ("mother: %g %g %g %g\n",vmother.Px(), vmother.Py(), vmother.Pz(), vmother.E()); + // printf ("vtot : %g %g %g %g\n",vtot.Px(), vtot.Py(), vtot.Pz(), vtot.E()); + // printf ("vvirt : %g %g %g %g\n",vvirt.Px(), vvirt.Py(), vvirt.Pz(), vvirt.E()); + // printf ("vdimu : %g %g %g %g\n",vdimu.Px(), vdimu.Py(), vdimu.Pz(), vdimu.E()); } @@ -616,7 +812,7 @@ Double_t AliGenMUONLMR::FormFactor(Double_t q2, Int_t decay){ Double_t lambda2inv = 0; switch (decay) { case 0: // eta -> mu mu gamma - // eta -> l+ l- gamma: pole approximation + // eta -> l+ l- gamma: pole approximation lambda2inv = 1.95; mass2 = fParticle[kEtaLMR]->GetMass() * fParticle[kEtaLMR]->GetMass(); if (q2 < mass2) ff2 = TMath::Power(1./(1.-lambda2inv*q2),2); @@ -639,7 +835,7 @@ Double_t AliGenMUONLMR::FormFactor(Double_t q2, Int_t decay){ else ff2 = 0; break; default: - printf ("FormFactor: Decay not found\n"); + AliError ("FormFactor: Decay not found"); return 0; break; } @@ -648,7 +844,7 @@ Double_t AliGenMUONLMR::FormFactor(Double_t q2, Int_t decay){ //____________________________________________________________ -Double_t AliGenMUONLMR::RhoLineShapeNew(const Double_t *x, const Double_t* /*para*/){ +Double_t AliGenMUONLMR::RhoLineShapeNew(Double_t *x, Double_t */*para*/){ //new parameterization implemented by Hiroyuki Sako (GSI) Double_t mass = *x; double r, GammaTot; @@ -658,6 +854,7 @@ Double_t AliGenMUONLMR::RhoLineShapeNew(const Double_t *x, const Double_t* /*par Double_t Gamma0 = TDatabasePDG::Instance()->GetParticle("rho0")->Width(); const double Norm = 0.0744416*1.01; + // 0.0744416 at m = 0.72297 // is the max number with Norm=1 (for rho) @@ -687,4 +884,3 @@ Double_t AliGenMUONLMR::RhoLineShapeNew(const Double_t *x, const Double_t* /*par return r; } -