X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=TDPMjet%2FAliGenDPMjet.cxx;h=4692ccf07907d226d265136157c061822d725050;hp=4e9cfc908ae07045919322f951355c2d5a6e0529;hb=d2593e7a88872ac14c838d07584053adb37839d0;hpb=1a52e0edb2bad3926b7e3a2784a35116d68abfc3 diff --git a/TDPMjet/AliGenDPMjet.cxx b/TDPMjet/AliGenDPMjet.cxx index 4e9cfc908ae..4692ccf0790 100644 --- a/TDPMjet/AliGenDPMjet.cxx +++ b/TDPMjet/AliGenDPMjet.cxx @@ -28,79 +28,112 @@ #include #include #include - +#include +#include "AliRunLoader.h" #include "AliGenDPMjet.h" #include "AliGenDPMjetEventHeader.h" -#include "AliPythia.h" #include "AliRun.h" +#include "AliDpmJetRndm.h" +#include "AliHeader.h" +#include "AliStack.h" +#include "AliMC.h" - - ClassImp(AliGenDPMjet) - +ClassImp(AliGenDPMjet) //______________________________________________________________________________ AliGenDPMjet::AliGenDPMjet() - :AliGenMC() + :AliGenMC(), + fBeamEn(2750.), + fEnergyCMS(5500.), + fMinImpactParam(0.), + fMaxImpactParam(5.), + fICentr(0), + fSelectAll(0), + fFlavor(0), + fTrials(0), + fSpectators(1), + fSpecn(0), + fSpecp(0), + fDPMjet(0), + fNoGammas(0), + fLHC(0), + fPi0Decay(0), + fGenImpPar(0.), + fProcess(kDpmMb) { // Constructor - fParticles = 0; - fDPMjet = 0; + AliDpmJetRndm::SetDpmJetRandom(GetRandom()); } //______________________________________________________________________________ AliGenDPMjet::AliGenDPMjet(Int_t npart) - :AliGenMC(npart) + :AliGenMC(npart), + fBeamEn(2750.), + fEnergyCMS(5500.), + fMinImpactParam(0.), + fMaxImpactParam(5.), + fICentr(0), + fSelectAll(0), + fFlavor(0), + fTrials(0), + fSpectators(1), + fSpecn(0), + fSpecp(0), + fDPMjet(0), + fNoGammas(0), + fLHC(0), + fPi0Decay(0), + fGenImpPar(0.), + fProcess(kDpmMb) { // Default PbPb collisions at 5. 5 TeV // fName = "DPMJET"; fTitle= "Particle Generator using DPMJET"; - - SetBeamEnergy(); - SetEnergyCMS(); SetTarget(); SetProjectile(); - SetCentral(); - SetImpactParameterRange(); - SetBoostLHC(); - - fKeep = 0; - fDecaysOff = 1; - fEvaluate = 0; - fSelectAll = 0; - fFlavor = 0; - fSpectators = 1; - fEventVertex.Set(3); - - fParticles = new TClonesArray("TParticle",10000); - - // Set random number generator - if (!sRandom) sRandom = fRandom; - fDPMjet = 0; - // Instance AliPythia - AliPythia::Instance(); - + fVertex.Set(3); + AliDpmJetRndm::SetDpmJetRandom(GetRandom()); } +AliGenDPMjet::AliGenDPMjet(const AliGenDPMjet &/*Dpmjet*/) + :AliGenMC(), + fBeamEn(2750.), + fEnergyCMS(5500.), + fMinImpactParam(0.), + fMaxImpactParam(5.), + fICentr(0), + fSelectAll(0), + fFlavor(0), + fTrials(0), + fSpectators(1), + fSpecn(0), + fSpecp(0), + fDPMjet(0), + fNoGammas(0), + fLHC(0), + fPi0Decay(0), + fGenImpPar(0.), + fProcess(kDpmMb) +{ + // Dummy copy constructor +} //______________________________________________________________________________ AliGenDPMjet::~AliGenDPMjet() { // Destructor - delete fParticles; } - - //______________________________________________________________________________ void AliGenDPMjet::Init() { // Initialization - SetMC(new TDPMjet(fAProjectile, fZProjectile, fATarget, fZTarget, + SetMC(new TDPMjet(fProcess, fAProjectile, fZProjectile, fATarget, fZTarget, fBeamEn,fEnergyCMS)); - fDPMjet=(TDPMjet*) fgMCEvGen; + fDPMjet=(TDPMjet*) fMCEvGen; // // **** Flag to force central production // fICentr=1. central production forced @@ -108,19 +141,15 @@ void AliGenDPMjet::Init() // fICentr<-99 -> fraction of x-sec. = XSFRAC // fICentr=-1. -> evaporation/fzc suppressed // fICentr<-1. -> evaporation/fzc suppressed - fDPMjet->SetfFCentr(fICentr); - fDPMjet->SetbRange(fMinImpactParam,fMaxImpactParam); + if (fAProjectile == 1 && fZProjectile == 1) fDPMjet->SetfIdp(1); + fDPMjet->SetfFCentr(fICentr); + fDPMjet->SetbRange(fMinImpactParam, fMaxImpactParam); + fDPMjet->SetPi0Decay(fPi0Decay); // // Initialize DPMjet // fDPMjet->Initialize(); - //if (fEvaluate) EvaluateCrossSections(); - - // Initialize AliIonPDGCodes to define ion PDG codes - /*AliIonPDGCodes *PDGcodes = new AliIonPDGCodes(); - PDGcodes->AddParticlesToPdgDataBase(); - PDGcodes->MapPDGGEant3Codes();*/ } @@ -131,31 +160,18 @@ void AliGenDPMjet::Generate() Float_t polar[3] = {0,0,0}; Float_t origin[3] = {0,0,0}; - Float_t origin0[3] = {0,0,0}; - Float_t p[3], random[6]; + Float_t p[3]; Float_t tof; // converts from mm/c to s const Float_t kconv = 0.001/2.999792458e8; Int_t nt = 0; Int_t jev = 0; - Int_t j, kf, ks, imo; + Int_t kf, ks, imo; kf = 0; - fTrials = 0; - for (j = 0;j < 3; j++) origin0[j] = fOrigin[j]; - if(fVertexSmear == kPerEvent) { - Float_t dv[3]; - dv[2] = 1.e10; - while(TMath::Abs(dv[2])>fCutVertexZ*fOsigma[2]) { - Rndm(random,6); - for (j=0; j < 3; j++) { - dv[j] = fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())* - TMath::Sqrt(-2*TMath::Log(random[2*j+1])); - } - } - for (j=0; j < 3; j++) origin0[j] += dv[j]; - } + // Set collision vertex position + if (fVertexSmear == kPerEvent) Vertex(); while(1) { @@ -167,13 +183,13 @@ void AliGenDPMjet::Generate() fDPMjet->GenerateEvent(); fTrials++; - fDPMjet->ImportParticles(fParticles,"All"); + fDPMjet->ImportParticles(&fParticles,"All"); if (fLHC) Boost(); // Temporaneo fGenImpPar = fDPMjet->GetBImpac(); - Int_t np = fParticles->GetEntriesFast(); + Int_t np = fParticles.GetEntriesFast(); printf("\n **************************************************%d\n",np); Int_t nc=0; if (np==0) continue; @@ -186,19 +202,10 @@ void AliGenDPMjet::Generate() pSelected[i] = 0; } -// Get event vertex - - //TParticle *iparticle = (TParticle*) fParticles->At(0); - fEventVertex[0] = origin0[0]; - fEventVertex[1] = origin0[1]; - fEventVertex[2] = origin0[2]; - //for(int jj=0; jj<3; jj++) printf(" fEventVertex[%d] = %f\n",jj,fEventVertex[jj]); - // First select parent particles for (i = 0; iAt(i); - //if(!iparticle) iparticle->Print(""); + TParticle *iparticle = (TParticle *) fParticles.At(i); // Is this a parent particle ? @@ -230,8 +237,7 @@ void AliGenDPMjet::Generate() for (i=0; iAt(i); - //if(!iparticle) printf(" i = %d -> iparticle==0\n",i); + TParticle *iparticle = (TParticle *) fParticles.At(i); // Is this a final state particle ? @@ -265,154 +271,50 @@ void AliGenDPMjet::Generate() // Write particles to stack for (i = 0; iAt(i); + TParticle * iparticle = (TParticle *) fParticles.At(i); Bool_t hasMother = (iparticle->GetFirstMother()>=0); - Bool_t hasDaughter = (iparticle->GetFirstDaughter()>=0); - if (pSelected[i]) { kf = iparticle->GetPdgCode(); - if(kf==99999) continue; ks = iparticle->GetStatusCode(); p[0] = iparticle->Px(); p[1] = iparticle->Py(); p[2] = iparticle->Pz(); - origin[0] = origin0[0]+iparticle->Vx()/10; - origin[1] = origin0[1]+iparticle->Vy()/10; - origin[2] = origin0[2]+iparticle->Vz()/10; + origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm] + origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm] + origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm] + tof = kconv*iparticle->T(); imo = -1; TParticle* mother = 0; if (hasMother) { imo = iparticle->GetFirstMother(); - mother = (TParticle *) fParticles->At(imo); - imo = (mother->GetPdgCode() != 92) ? imo = newPos[imo] : -1; + mother = (TParticle *) fParticles.At(imo); + imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1; } // if has mother - // --- Nucleons and photons from nucleus evaporation - //if(ks==-1 && (kf==2212 || kf==2112 || kf==22)) imo=-1; - - // --- Offset added to nuclei, alpha particles, deuteron - // --- and triton PDG codes (according to TGeant3 PDG codes) - //if(ks==1000 || ks==1001) kf += 10000000; - if(kf>10000) kf += 10000000; - - if(kf>10000 || ks==1000 || ks==1001) printf(" kf = %d, ks = %d imo = %d \n", kf,ks,imo); - - Bool_t tFlag = (fTrackIt && (ks == 1)); -/* // --- Particle NOT to be tracked: - // --- (1) chains from valence quark system (idhkk=99999) - // --- (2) quark, diquarks and gluons (idhkk=1->8,...,idhkk=21) - if(kf==99999 || kf==1 || kf==2 || kf==3 || kf==4 || kf==5 || - kf==6 || kf==7 || kf==8 || kf==21 || - kf==1103 || kf==2101 || kf==2103 || kf==2203 || kf==3101 || - kf==3103 || kf==3201 || kf==3203 || kf==3303 || kf==4101 || - kf==4103 || kf==4201 || kf==4203 || kf==4301 || kf==4303 || - kf==4403 || kf==5101 || kf==5103 || kf==5201 || kf==5203 || - kf==5301 || kf==5303 || kf==5401 || kf==5403 || kf==5503) - tFlag=kFALSE; -*/ - SetTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks); + PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks); KeepTrack(nt); newPos[i] = nt; } // if selected } // particle loop delete[] newPos; delete[] pSelected; - - printf("\n I've put %i particles on the stack \n",nc); if (nc>0) { jev += nc; if (jev >= fNpart || fNpart == -1) { - printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev); break; } - printf("\n Not breaking ### Trials: %i %i %i\n",fTrials, fNpart, jev); } } // event loop MakeHeader(); SetHighWaterMark(nt); } - -//______________________________________________________________________________ -void AliGenDPMjet::KeepFullEvent() -{ - fKeep=1; -} - - -//______________________________________________________________________________ -/*void AliGenDPMjet::EvaluateCrossSections() -{ -// Glauber Calculation of geometrical x-section -// - Float_t xTot = 0.; // barn - Float_t xTotHard = 0.; // barn - Float_t xPart = 0.; // barn - Float_t xPartHard = 0.; // barn - Float_t sigmaHard = 0.1; // mbarn - Float_t bMin = 0.; - Float_t bMax = fDPMjet->GetProjRadius()+fDPMjet->GetTargRadius(); - const Float_t kdib = 0.2; - Int_t kMax = Int_t((bMax-bMin)/kdib)+1; - - - printf("\n Projectile Radius (fm): %f \n",fDPMjet->GetProjRadius()); - printf("\n Target Radius (fm): %f \n",fDPMjet->GetTargRadius()); - Int_t i; - Float_t oldvalue= 0.; - - Float_t* b = new Float_t[kMax]; - Float_t* si1 = new Float_t[kMax]; - Float_t* si2 = new Float_t[kMax]; - - for (i = 0; i < kMax; i++) - { - Float_t xb = bMin+i*kdib; - Float_t ov; - ov=fDPMjet->Profile(xb); - // ATT!->Manca la x-sec anel. nucleone-nucleone - Float_t gb = 2.*0.01*fDPMjet->TMath::Pi()*kdib*xb*(1.-TMath::Exp(-fDPMjet->GetXSFrac()*ov)); - Float_t gbh = 2.*0.01*fDPMjet->TMath::Pi()*kdib*xb*sigmaHard*ov; - xTot+=gb; - xTotHard += gbh; - if (xb > fMinImpactParam && xb < fMaxImpactParam) - { - xPart += gb; - xPartHard += gbh; - } - - if(oldvalue) if ((xTot-oldvalue)/oldvalue<0.0001) break; - oldvalue = xTot; - printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot); - printf("\n Hard cross section (barn): %d %f %f \n\n",i, xb, xTotHard); - if (i>0) { - si1[i] = gb/kdib; - si2[i] = gbh/gb; - b[i] = xb; - } - } - - printf("\n Total cross section (barn): %f \n",xTot); - printf("\n Hard cross section (barn): %f \n \n",xTotHard); - printf("\n Partial cross section (barn): %f %f \n",xPart, xPart/xTot*100.); - printf("\n Partial hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.); - -// Store result as a graph - b[0] = 0; - si1[0] = 0; - si2[0]=si2[1]; - - fDsigmaDb = new TGraph(i, b, si1); - fDnDb = new TGraph(i, b, si2); -}*/ - - //______________________________________________________________________________ Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle) { @@ -429,7 +331,7 @@ Bool_t AliGenDPMjet::DaughtersSelection(TParticle* iparticle) imin = iparticle->GetFirstDaughter(); imax = iparticle->GetLastDaughter(); for (i = imin; i <= imax; i++){ - TParticle * jparticle = (TParticle *) fParticles->At(i); + TParticle * jparticle = (TParticle *) fParticles.At(i); Int_t ip = jparticle->GetPdgCode(); if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) { selected=kTRUE; break; @@ -480,42 +382,6 @@ Bool_t AliGenDPMjet::Stable(TParticle* particle) } - - -//______________________________________________________________________________ -/*void AliGenDPMjet::Boost(TClonesArray* particles) -{ -// -// Boost cms into LHC lab frame -// - Double_t dy = - 0.5 * TMath::Log(Double_t(fZProjectile) * Double_t(fATarget) / - (Double_t(fZTarget) * Double_t(fAProjectile))); - Double_t beta = TMath::TanH(dy); - Double_t gamma = 1./TMath::Sqrt(1.-beta*beta); - Double_t gb = gamma * beta; - - printf("\n Boosting particles to lab frame %f %f %f", dy, beta, gamma); - - Int_t i; - Int_t np = particles->GetEntriesFast(); - for (i = 0; i < np; i++) - { - TParticle* iparticle = (TParticle*) fParticles->At(i); - - Double_t e = iparticle->Energy(); - Double_t px = iparticle->Px(); - Double_t py = iparticle->Py(); - Double_t pz = iparticle->Pz(); - - Double_t eb = gamma * e - gb * pz; - Double_t pzb = -gb * e + gamma * pz; - - iparticle->SetMomentum(px, py, pzb, eb); - } -}*/ - - - //______________________________________________________________________________ void AliGenDPMjet::MakeHeader() { @@ -524,65 +390,42 @@ void AliGenDPMjet::MakeHeader() ((AliGenDPMjetEventHeader*) header)->SetNProduced(fDPMjet->GetNumStablePc()); ((AliGenDPMjetEventHeader*) header)->SetImpactParameter(fDPMjet->GetBImpac()); ((AliGenDPMjetEventHeader*) header)->SetTotalEnergy(fDPMjet->GetTotEnergy()); -// ((AliGenDPMjetEventHeader*) header)->SetWounded(fDPMjet->GetProjWounded(), -// fDPMjet->GetTargWounded()); ((AliGenDPMjetEventHeader*) header)->SetParticipants(fDPMjet->GetfIp(), fDPMjet->GetfIt()); - /*((AliGenDPMjetEventHeader*) header)->SetCollisions(fDPMjet->GetN0(), - fDPMjet->GetN01(), - fDPMjet->GetN10(), - fDPMjet->GetN11());*/ - ((AliGenDPMjetEventHeader*) header)->SetSpectators(fSpecn, fSpecp); - + ((AliGenDPMjetEventHeader*) header)->SetProcessType(fDPMjet->GetProcessCode()); // Bookkeeping for kinematic bias ((AliGenDPMjetEventHeader*) header)->SetTrials(fTrials); // Event Vertex - header->SetPrimaryVertex(fEventVertex); + header->SetPrimaryVertex(fVertex); gAlice->SetGenEventHeader(header); + AddHeader(header); } +void AliGenDPMjet::AddHeader(AliGenEventHeader* header) +{ + // Add header to container or runloader + if (fContainer) { + fContainer->AddHeader(header); + } else { + AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header); + } +} //______________________________________________________________________________ -AliGenDPMjet& AliGenDPMjet::operator=(const AliGenDPMjet& rhs) +AliGenDPMjet& AliGenDPMjet::operator=(const AliGenDPMjet& /*rhs*/) { // Assignment operator return *this; } - -//______________________________________________________________________________ -#ifndef WIN32 -# define dt_rndm_dpmjet dt_rndm_dpmjet_ -# define dt_rndmst_dpmjet dt_rndmst_dpmjet_ -# define dt_rndmin_dpmjet dt_rndmin_dpmjet_ -# define dt_rndmou_dpmjet dt_rndmou_dpmjet_ -# define type_of_call -#else -# define dt_rndm_dpmjet DT_RNDM_DPMJET -# define dt_rndmst_dpmjet DT_RNDMST_DPMJET -# define dt_rndmin_dpmjet DT_RNDMIN_DPMJET -# define dt_rndmou_dpmjet DT_RNDMOU_DPMJET -# define type_of_call _stdcall -#endif - - -extern "C" { - void type_of_call dt_rndmst_dpmjet(Int_t &, Int_t &, Int_t &, Int_t &) - {printf("Dummy version of dt_rndmst reached\n");} - - void type_of_call dt_rndmin_dpmjet(Int_t &, Int_t &, Int_t &, Int_t &, Int_t &, Int_t &) - {printf("Dummy version of dt_rndmin reached\n");} - - void type_of_call dt_rndmou_dpmjet(Int_t &, Int_t &, Int_t &, Int_t &, Int_t &, Int_t &) - {printf("Dummy version of dt_rndmou reached\n");} - - Double_t type_of_call dt_rndm_dpmjet(Int_t &) - { - Float_t r; - do r = sRandom->Rndm(); while(0 >= r || r >= 1); - return r; - } +void AliGenDPMjet::FinishRun() +{ + // Print run statistics + fDPMjet->Dt_Dtuout(); } + + +//______________________________________________________________________________