From: morsch Date: Fri, 13 Jul 2001 10:58:54 +0000 (+0000) Subject: - Some coded moved to AliGenMC X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=fff02fee2c2b64f6690f1cba68f5bb819951a093;p=u%2Fmrichter%2FAliRoot.git - Some coded moved to AliGenMC - Improved handling of secondary vertices. --- diff --git a/EVGEN/AliGenParam.cxx b/EVGEN/AliGenParam.cxx index 80dc91c5e01..07dcca9c044 100644 --- a/EVGEN/AliGenParam.cxx +++ b/EVGEN/AliGenParam.cxx @@ -15,6 +15,9 @@ /* $Log$ +Revision 1.30 2001/06/15 07:55:04 morsch +Put only first generation decay products on the stack. + Revision 1.29 2001/03/27 10:58:41 morsch Initialize decayer before generation. Important if run inside cocktail. @@ -111,18 +114,13 @@ AliGenParam::AliGenParam() fYPara = 0; fParam = 0; fAnalog = kAnalog; - SetCutOnChild(); - SetChildMomentumRange(); - SetChildPtRange(); - SetChildPhiRange(); - SetChildThetaRange(); SetDeltaPt(); // // Set random number generator sRandom = fRandom; } -AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library, Int_t param, char* tname):AliGenerator(npart) +AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library, Int_t param, char* tname):AliGenMC(npart) { // Constructor using number of particles parameterisation id and library @@ -134,14 +132,7 @@ AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library, Int_t param, char* t fYPara = 0; fParam = param; fAnalog = kAnalog; - fChildSelect.Set(5); - for (Int_t i=0; i<5; i++) fChildSelect[i]=0; SetForceDecay(); - SetCutOnChild(); - SetChildMomentumRange(); - SetChildPtRange(); - SetChildPhiRange(); - SetChildThetaRange(); SetDeltaPt(); // // Set random number generator @@ -150,7 +141,7 @@ AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library, Int_t param, char* t //____________________________________________________________ -AliGenParam::AliGenParam(Int_t npart, Int_t param, char* tname):AliGenerator(npart) +AliGenParam::AliGenParam(Int_t npart, Int_t param, char* tname):AliGenMC(npart) { // Constructor using parameterisation id and number of particles // @@ -179,7 +170,7 @@ AliGenParam::AliGenParam(Int_t npart, Int_t param, Double_t (*PtPara) (Double_t*, Double_t*), Double_t (*YPara ) (Double_t* ,Double_t*), Int_t (*IpPara) (TRandom *)) - :AliGenerator(npart) + :AliGenMC(npart) { // Constructor // Gines Martinez 1/10/99 @@ -265,34 +256,7 @@ void AliGenParam::Init() fDecayer->Init(); // - switch (fForceDecay) - { - case kSemiElectronic: - case kDiElectron: - case kBJpsiDiElectron: - case kBPsiPrimeDiElectron: - fChildSelect[0]=11; - break; - case kSemiMuonic: - case kDiMuon: - case kBJpsiDiMuon: - case kBPsiPrimeDiMuon: - fChildSelect[0]=13; - break; - case kPiToMu: - fChildSelect[0]=13; - break; - case kKaToMu: - fChildSelect[0]=13; - break; - case kHadronicD: -// Implement me !! - break; - case kNoDecay: - break; - case kAll: - break; - } + AliGenMC::Init(); } //____________________________________________________________ @@ -317,14 +281,14 @@ void AliGenParam::Generate() Float_t pt, pl, ptot; // Transverse, logitudinal and total momenta of the parent particle Float_t phi, theta; // Phi and theta spherical angles of the parent particle momentum Float_t p[3], pc[3], - och[3], pch[10][3]; // Momentum, polarisation and origin of the children particles from lujet + och[3]; // Momentum, polarisation and origin of the children particles from lujet Float_t ty, xmt; - Int_t nt, i, j, kfch[10]; + Int_t nt, i, j; Float_t wgtp, wgtch; Double_t dummy; static TClonesArray *particles; // - if(!particles) particles=new TClonesArray("TParticle",1000); + if(!particles) particles = new TClonesArray("TParticle",1000); static TDatabasePDG *pDataBase = new TDatabasePDG(); if(!pDataBase) pDataBase = new TDatabasePDG(); @@ -334,18 +298,12 @@ void AliGenParam::Generate() // Calculating vertex position per event for (j=0;j<3;j++) origin0[j]=fOrigin[j]; if(fVertexSmear==kPerEvent) { -// Rndm(random,6); -// for (j=0;j<3;j++) { -// origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())* -// TMath::Sqrt(-2*TMath::Log(random[2*j+1])); -// TMath::Sqrt(-2*TMath::Log(random[2*j+1])); -// } -// } Vertex(); for (j=0;j<3;j++) origin0[j]=fVertex[j]; } Int_t ipa=0; + // Generating fNpart particles while (ipafPMax) continue; +// p[0]=pt*TMath::Cos(phi); p[1]=pt*TMath::Sin(phi); p[2]=pl; @@ -403,6 +362,7 @@ void AliGenParam::Generate() // if fForceDecay == none Primary particle is tracked by GEANT // (In the latest, make sure that GEANT actually does all the decays you want) // + if (fForceDecay != kNoDecay) { // Using lujet to decay particle Float_t energy=TMath::Sqrt(ptot*ptot+am*am); @@ -412,50 +372,71 @@ void AliGenParam::Generate() // select decay particles Int_t np=fDecayer->ImportParticles(particles); Int_t ncsel=0; - + Int_t* pFlag = new Int_t[np]; + Int_t* pParent = new Int_t[np]; + Int_t* pSelected = new Int_t[np]; + Int_t* trackIt = new Int_t[np]; + + for (i=0; i1) { TParticle* iparticle = (TParticle *) particles->At(0); - Int_t ipF = iparticle->GetFirstDaughter(); - Int_t ipL = iparticle->GetLastDaughter(); - for (i = ipF-1; iAt(i); Int_t kf = iparticle->GetPdgCode(); + Int_t ks = iparticle->GetStatusCode(); +// flagged particle + + if (pFlag[i] == 1) { + printf("\n deselected %d", i); + ipF = iparticle->GetFirstDaughter(); + ipL = iparticle->GetLastDaughter(); + if (ipF > 0) for (j=ipF-1; j .3 mum) + + if (ks != 1) { + TParticlePDG *particle = pDataBase->GetParticle(kf); + + Double_t lifeTime = fDecayer->GetLifetime(kf); + Double_t mass = particle->Mass(); + Double_t width = particle->Width(); + printf("\n lifetime %d %e %e %e ", i, lifeTime, mass, width); + if (lifeTime > (Double_t) 1.e-15) { + ipF = iparticle->GetFirstDaughter(); + ipL = iparticle->GetLastDaughter(); + if (ipF > 0) for (j=ipF-1; jPx(); - pc[1]=iparticle->Py(); - pc[2]=iparticle->Pz(); - och[0]=origin0[0]+iparticle->Vx()/10; - och[1]=origin0[1]+iparticle->Vy()/10; - och[2]=origin0[2]+iparticle->Vz()/10; if (fCutOnChild) { - Float_t ptChild=TMath::Sqrt(pc[0]*pc[0]+pc[1]*pc[1]); - Float_t pChild=TMath::Sqrt(ptChild*ptChild+pc[2]*pc[2]); - Float_t thetaChild=TMath::ATan2(ptChild,pc[2]); - Float_t phiChild=TMath::ATan2(pc[1],pc[0]); - Bool_t childok = - ((ptChild > fChildPtMin && ptChild fChildPMin && pChild fChildThetaMin && thetaChild fChildPhiMin && phiChild Px(); + pc[1]=iparticle->Py(); + pc[2]=iparticle->Pz(); + Bool_t childok = KinematicSelection(iparticle, 1); + if(childok) { + pSelected[i] = 1; ncsel++; } else { ncsel=-1; break; } // child kine cuts } else { - pch[ncsel][0]=pc[0]; - pch[ncsel][1]=pc[1]; - pch[ncsel][2]=pc[2]; - kfch[ncsel]=kf; + pSelected[i] = 1; ncsel++; } // if child selection } // select muon @@ -466,18 +447,45 @@ void AliGenParam::Generate() if ((fCutOnChild && ncsel >0) || !fCutOnChild){ ipa++; // -// parent - gAlice-> - SetTrack(0,-1,iPart,p,origin0,polar,0,kPPrimary,nt,wgtp); - iparent=nt; +// Parent + gAlice->SetTrack(0, -1, iPart, p, origin0, polar, 0, kPPrimary, nt, wgtp); + pParent[0] = nt; gAlice->KeepTrack(nt); - for (i=0; i< ncsel; i++) { - gAlice->SetTrack(fTrackIt,iparent,kfch[i], - &pch[i][0],och,polar, - 0,kPDecay,nt,wgtch); - gAlice->KeepTrack(nt); +// +// Decay Products +// + for (i = 1; i < np; i++) { + if (pSelected[i]) { + TParticle* iparticle = (TParticle *) particles->At(i); + Int_t kf = iparticle->GetPdgCode(); + Int_t ipa = iparticle->GetFirstMother()-1; + + och[0] = origin0[0]+iparticle->Vx()/10; + och[1] = origin0[1]+iparticle->Vy()/10; + och[2] = origin0[2]+iparticle->Vz()/10; + pc[0] = iparticle->Px(); + pc[1] = iparticle->Py(); + pc[2] = iparticle->Pz(); + + if (ipa > -1) { + iparent = pParent[ipa]; + } else { + iparent = -1; + } + printf("\n SetTrack %d %d %d %d", i, kf, iparent, trackIt[i]); + gAlice->SetTrack(fTrackIt*trackIt[i], iparent, kf, + pc, och, polar, + 0, kPDecay, nt, wgtch); + pParent[i] = nt; + gAlice->KeepTrack(nt); + } } } // Decays by Lujet + + if (pFlag) delete[] pFlag; + if (pParent) delete[] pParent; + if (pSelected) delete[] pSelected; + if (trackIt) delete[] trackIt; } // kinematic selection else // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons) { @@ -489,48 +497,8 @@ void AliGenParam::Generate() } // while } // event loop gAlice->SetHighWaterMark(nt); - } -Bool_t AliGenParam::ChildSelected(Int_t ip) -{ -// True if particle is in list of selected children - for (Int_t i=0; i<5; i++) - { - if (fChildSelect[i]==ip) return kTRUE; - } - return kFALSE; -} - -Bool_t AliGenParam::KinematicSelection(TParticle *particle) -{ -// Perform kinematic cuts - Float_t px=particle->Px(); - Float_t py=particle->Py(); - Float_t pz=particle->Pz(); -// -// momentum cut - Float_t p=TMath::Sqrt(px*px+py*py+pz*pz); - if (p > fPMax || p < fPMin) - { -// printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax); - return kFALSE; - } - Float_t pt=TMath::Sqrt(px*px+py*py); - -// -// theta cut - Float_t theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(p))); - if (theta > fThetaMax || theta < fThetaMin) - { -// printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax); - return kFALSE; - } - - return kTRUE; -} - - AliGenParam& AliGenParam::operator=(const AliGenParam& rhs) { // Assignment operator diff --git a/EVGEN/AliGenParam.h b/EVGEN/AliGenParam.h index ec57eec89e9..dfd4dc8ef68 100644 --- a/EVGEN/AliGenParam.h +++ b/EVGEN/AliGenParam.h @@ -5,8 +5,7 @@ /* $Id$ */ -#include "AliGenerator.h" -#include "AliDecayer.h" +#include "AliGenMC.h" #include class AliPythia; @@ -16,7 +15,7 @@ class TF1; typedef enum { kAnalog, kNonAnalog} Weighting_t; //------------------------------------------------------------- -class AliGenParam : public AliGenerator +class AliGenParam : public AliGenMC { public: AliGenParam(); @@ -34,19 +33,7 @@ class AliGenParam : public AliGenerator // select particle type virtual void SetParam(Int_t param) {fParam = param;} // force decay type - virtual void SetForceDecay(Decay_t decay = kDiMuon) {fForceDecay = decay;} virtual void SetWeighting(Weighting_t flag = kAnalog) {fAnalog = flag;} - virtual void SetCutOnChild(Int_t flag = 0) {fCutOnChild = flag;} - virtual void SetChildMomentumRange(Float_t pmin = 0, Float_t pmax = 1.e10) - {fChildPMin = pmin; fChildPMax = pmax;} - virtual void SetChildPtRange(Float_t ptmin = 0, Float_t ptmax = 20.) - {fChildPtMin = ptmin; fChildPtMax = ptmax;} - virtual void SetChildPhiRange(Float_t phimin = -180., Float_t phimax = 180) - {fChildPhiMin = TMath::Pi()*phimin/180; - fChildPhiMax = TMath::Pi()*phimax/180;} - virtual void SetChildThetaRange(Float_t thetamin = 0, Float_t thetamax = 180) - {fChildThetaMin = TMath::Pi()*thetamin/180; - fChildThetaMax = TMath::Pi()*thetamax/180;} virtual void SetDeltaPt(Float_t delta=0.01) {fDeltaPt = delta;} AliGenParam & operator=(const AliGenParam & rhs); @@ -62,25 +49,8 @@ class AliGenParam : public AliGenerator Float_t fPtWgt; // Pt-weight Float_t fBias; // Biasing factor Int_t fTrials; // Number of trials - Decay_t fForceDecay; // Decay channel forced - Int_t fCutOnChild; // Cuts on decay products (children) are enabled/disabled - Float_t fChildPtMin; // Children minimum pT - Float_t fChildPtMax; // Children maximum pT - Float_t fChildPMin; // Children minimum p - Float_t fChildPMax; // Children maximum p - Float_t fChildPhiMin; // Children minimum phi - Float_t fChildPhiMax; // Children maximum phi - Float_t fChildThetaMin;// Children minimum theta - Float_t fChildThetaMax;// Children maximum theta Float_t fDeltaPt; // pT sampling in steps of fDeltaPt - TArrayI fChildSelect; // Children to be selected from decay products AliDecayer *fDecayer; // ! Pointer to pythia object for decays - private: - // check if particle is selected as child - Bool_t ChildSelected(Int_t ip); - // all kinematic selection goes here - Bool_t KinematicSelection(TParticle *particle); - ClassDef(AliGenParam,1) // Generator using parameterised pt- and y-distribution }; #endif diff --git a/EVGEN/AliGenPythia.cxx b/EVGEN/AliGenPythia.cxx index 418910ff549..d4cf4fc270a 100644 --- a/EVGEN/AliGenPythia.cxx +++ b/EVGEN/AliGenPythia.cxx @@ -15,6 +15,10 @@ /* $Log$ +Revision 1.37 2001/06/28 11:17:28 morsch +SetEventListRange setter added. Events in specified range are listed for +debugging. (Yuri Kharlov) + Revision 1.36 2001/03/30 07:05:49 morsch Final print-out in finish run. Write parton system for jet-production (preliminary solution). @@ -102,6 +106,7 @@ Introduction of the Copyright and cvs Log */ #include "AliGenPythia.h" +#include "AliGenPythiaEventHeader.h" #include "AliDecayerPythia.h" #include "AliRun.h" #include "AliPythia.h" @@ -112,7 +117,7 @@ Introduction of the Copyright and cvs Log ClassImp(AliGenPythia) AliGenPythia::AliGenPythia() - :AliGenerator() + :AliGenMC() { // Default Constructor fDecayer = new AliDecayerPythia(); @@ -120,7 +125,7 @@ AliGenPythia::AliGenPythia() } AliGenPythia::AliGenPythia(Int_t npart) - :AliGenerator(npart) + :AliGenMC(npart) { // default charm production at 5. 5 TeV // semimuonic decay @@ -129,9 +134,6 @@ AliGenPythia::AliGenPythia(Int_t npart) fXsection = 0.; fNucA1=0; fNucA2=0; - fParentSelect.Set(5); - fChildSelect.Set(5); - for (Int_t i=0; i<5; i++) fParentSelect[i]=fChildSelect[i]=0; SetProcess(); SetStrucFunc(); SetForceDecay(); @@ -141,6 +143,9 @@ AliGenPythia::AliGenPythia(Int_t npart) // Set random number generator sRandom=fRandom; SetEventListRange(); + fFlavorSelect = 0; + // Produced particles + fParticles = new TClonesArray("TParticle",1000); } AliGenPythia::AliGenPythia(const AliGenPythia & Pythia) @@ -171,10 +176,6 @@ void AliGenPythia::Init() fParentWeight=1./Float_t(fNpart); // // Forward Paramters to the AliPythia object - // gSystem->Exec("ln -s $ALICE_ROOT/data/Decay.table fort.1"); - // fPythia->Pyupda(2,1); - // gSystem->Exec("rm fort.1"); - fDecayer->SetForceDecay(fForceDecay); fDecayer->Init(); @@ -190,63 +191,49 @@ void AliGenPythia::Init() switch (fProcess) { case kPyCharm: - fParentSelect[0]=411; - fParentSelect[1]=421; - fParentSelect[2]=431; - fParentSelect[3]=4122; + fParentSelect[0] = 411; + fParentSelect[1] = 421; + fParentSelect[2] = 431; + fParentSelect[3] = 4122; + fFlavorSelect = 4; break; case kPyCharmUnforced: - fParentSelect[0]=411; - fParentSelect[1]=421; - fParentSelect[2]=431; - fParentSelect[3]=4122; + fParentSelect[0] = 411; + fParentSelect[1] = 421; + fParentSelect[2] = 431; + fParentSelect[3]= 4122; + fFlavorSelect = 4; break; case kPyBeauty: - fParentSelect[0]=511; - fParentSelect[1]=521; - fParentSelect[2]=531; - fParentSelect[3]=5122; + fParentSelect[0]= 511; + fParentSelect[1]= 521; + fParentSelect[2]= 531; + fParentSelect[3]= 5122; + fParentSelect[4]= 5132; + fParentSelect[5]= 5232; + fParentSelect[6]= 5332; + fFlavorSelect = 5; break; case kPyBeautyUnforced: - fParentSelect[0]=511; - fParentSelect[1]=521; - fParentSelect[2]=531; - fParentSelect[3]=5122; + fParentSelect[0] = 511; + fParentSelect[1] = 521; + fParentSelect[2] = 531; + fParentSelect[3] = 5122; + fParentSelect[4] = 5132; + fParentSelect[5] = 5232; + fParentSelect[6] = 5332; + fFlavorSelect = 5; break; case kPyJpsiChi: case kPyJpsi: - fParentSelect[0]=443; + fParentSelect[0] = 443; break; case kPyMb: case kPyJets: case kPyDirectGamma: break; } - - switch (fForceDecay) - { - case kSemiElectronic: - case kDiElectron: - case kBJpsiDiElectron: - case kBPsiPrimeDiElectron: - fChildSelect[0]=kElectron; - break; - case kSemiMuonic: - case kDiMuon: - case kBJpsiDiMuon: - case kBPsiPrimeDiMuon: - case kPiToMu: - case kKaToMu: - fChildSelect[0]=kMuonMinus; - break; - case kHadronicD: - fChildSelect[0]=kPiPlus; - fChildSelect[1]=kKPlus; - break; - case kAll: - case kNoDecay: - break; - } + AliGenMC::Init(); } void AliGenPythia::Generate() @@ -256,138 +243,147 @@ void AliGenPythia::Generate() Float_t polar[3] = {0,0,0}; Float_t origin[3] = {0,0,0}; - Float_t originP[3] = {0,0,0}; Float_t origin0[3] = {0,0,0}; - Float_t p[3], pP[4]; -// Float_t random[6]; - static TClonesArray *particles; + Float_t p[3]; // converts from mm/c to s const Float_t kconv=0.001/2.999792458e8; - - // Int_t nt=0; - Int_t ntP=0; Int_t jev=0; Int_t j, kf; - - if(!particles) particles=new TClonesArray("TParticle",1000); - fTrials=0; + +// Set collision vertex position for (j=0;j<3;j++) origin0[j]=fOrigin[j]; if(fVertexSmear==kPerEvent) { fPythia->SetMSTP(151,1); for (j=0;j<3;j++) { - fPythia->SetPARP(151+j, fOsigma[j]/10.); + fPythia->SetPARP(151+j, fOsigma[j]*10.); } } else if (fVertexSmear==kPerTrack) { fPythia->SetMSTP(151,0); } - +// event loop while(1) { fPythia->Pyevnt(); if (gAlice->GetEvNumber()>=fDebugEventFirst && gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1); fTrials++; - fPythia->ImportParticles(particles,"All"); - Int_t np = particles->GetEntriesFast(); + fPythia->ImportParticles(fParticles,"All"); + +// +// +// + Int_t i; + + Int_t np = fParticles->GetEntriesFast(); + Int_t* pParent = new Int_t[np]; + Int_t* pSelected = new Int_t[np]; + Int_t* trackIt = new Int_t[np]; + for (i=0; i< np-1; i++) { + pParent[i] = -1; + pSelected[i] = 0; + } printf("\n **************************************************%d\n",np); - Int_t nc=0; + Int_t nc = 0; if (np == 0 ) continue; if (fProcess != kPyMb && fProcess != kPyJets && fProcess != kPyDirectGamma) { - for (Int_t i = 0; iAt(i); + for (i = 0; iAt(i); Int_t ks = iparticle->GetStatusCode(); kf = CheckPDGCode(iparticle->GetPdgCode()); +// No initial state partons if (ks==21) continue; - fChildWeight=(fDecayer->GetPartialBranchingRatio(kf))*fParentWeight; // -// Parent - if (ParentSelected(TMath::Abs(kf))) { - if (KinematicSelection(iparticle)) { - if (nc==0) { +// Heavy Flavor Selection // -// Store information concerning the hard scattering process + // quark ? + kf = TMath::Abs(kf); + Int_t kfl = kf; + // meson ? + if (kfl > 10) kfl/=100; + // baryon + if (kfl > 10) kfl/=10; + if (kfl > 10) kfl/=10; + + Int_t ipa = iparticle->GetFirstMother()-1; + Int_t kfMo = 0; + + if (ipa > -1) { + TParticle * mother = (TParticle *) fParticles->At(ipa); + kfMo = TMath::Abs(mother->GetPdgCode()); + } +// printf("\n particle (all) %d %d %d", i, pSelected[i], kf); + if (kfl >= fFlavorSelect) { // - Float_t massP = fPythia->GetPARI(13); - Float_t ptP = fPythia->GetPARI(17); - Float_t yP = fPythia->GetPARI(37); - Float_t xmtP = sqrt(ptP*ptP+massP*massP); - Float_t ty = Float_t(TMath::TanH(yP)); - pP[0] = ptP; - pP[1] = 0; - pP[2] = xmtP*ty/sqrt(1.-ty*ty); - pP[3] = massP; - gAlice->SetTrack(0,-1,-1, - pP,originP,polar, - 0,kPPrimary,ntP,fParentWeight); -// 0,"Hard Scat.",ntP,fParentWeight); - gAlice->KeepTrack(ntP); - } - nc++; +// Heavy flavor hadron or quark // -// store parent track information - 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.; - - Int_t ifch=iparticle->GetFirstDaughter(); - Int_t ilch=iparticle->GetLastDaughter(); - - if ((ifch !=0 && ilch !=0) || fForceDecay == kNoDecay) { - Int_t trackit=0; - if (fForceDecay == kNoDecay) trackit = 1; - gAlice->SetTrack(trackit,ntP,kf, - p,origin,polar, - 0,kPPrimary,nt,fParentWeight); - gAlice->KeepTrack(nt); - Int_t iparent = nt; +// Kinematic seletion on final state heavy flavor mesons + if (ParentSelected(kf) && !KinematicSelection(iparticle, 0)) + { + nc = -1; + break; + } + pSelected[i] = 1; +// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf); + } else { +// Kinematic seletion on decay products + if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf) + && !KinematicSelection(iparticle, 1)) + { + nc = -1; + break; + } // -// Children - if (fForceDecay != kNoDecay) { - for (j=ifch; j<=ilch; j++) - { - TParticle * ichild = - (TParticle *) particles->At(j-1); - kf = CheckPDGCode(ichild->GetPdgCode()); +// Decay products +// Select if mother was selected and is not tracked + + if (pSelected[ipa] && + !trackIt[ipa] && // mother will be tracked ? + kfMo != 5 && // mother is b-quark, don't store fragments + kfMo != 4 && // mother is c-quark, don't store fragments + kf != 92) // don't store string + { // +// Semi-stable or de-selected: diselect decay products: // - if (ChildSelected(TMath::Abs(kf))) { - origin[0]=origin0[0]+ichild->Vx()/10.; - origin[1]=origin0[1]+ichild->Vy()/10.; - origin[2]=origin0[2]+ichild->Vz()/10.; - p[0]=ichild->Px(); - p[1]=ichild->Py(); - p[2]=ichild->Pz(); - Float_t tof=kconv*ichild->T(); - gAlice->SetTrack(fTrackIt, iparent, kf, - p,origin,polar, - tof,kPDecay,nt,fChildWeight); - gAlice->KeepTrack(nt); - } // select child - } // child loop - } - } - } // kinematic selection - } // select particle - } // particle loop - } else { - for (Int_t i = 0; iAt(i); - kf = CheckPDGCode(iparticle->GetPdgCode()); - Int_t ks = iparticle->GetStatusCode(); - Int_t km = iparticle->GetFirstMother(); - // printf("\n process %d %d\n", ks,km); - - if ((ks==1 && kf!=0 && KinematicSelection(iparticle)) || - (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) { - nc++; // -// store track information + if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > 10e-15) + { + Int_t ipF = iparticle->GetFirstDaughter(); + Int_t ipL = iparticle->GetLastDaughter(); + if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1; + } +// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf); + pSelected[i] = (pSelected[i] == -1) ? 0 : 1; + } + } + if (pSelected[i] == -1) pSelected[i] = 0; + if (!pSelected[i]) continue; + nc++; +// Decision on tracking + trackIt[i] = 0; +// +// Track final state particle + if (ks == 1) trackIt[i] = 1; +// Track semi-stable particles + if ((ks ==1) || (fDecayer->GetLifetime(kf)> 10e-15)) trackIt[i] = 1; +// Track particles selected by process if undecayed. + if (fForceDecay == kNoDecay) { + if (ParentSelected(kf)) trackIt[i] = 1; + } else { + if (ParentSelected(kf)) trackIt[i] = 0; + } +// +// + + } // particle selection loop + if (nc > -1) { + for (i = 0; iAt(i); + kf = CheckPDGCode(iparticle->GetPdgCode()); p[0]=iparticle->Px(); p[1]=iparticle->Py(); p[2]=iparticle->Pz(); @@ -395,13 +391,19 @@ void AliGenPythia::Generate() origin[1]=origin0[1]+iparticle->Vy()/10.; origin[2]=origin0[2]+iparticle->Vz()/10.; Float_t tof=kconv*iparticle->T(); - gAlice->SetTrack(fTrackIt,-1,kf,p,origin,polar, - tof,kPPrimary,nt); - gAlice->KeepTrack(nt); - } // select particle - } // particle loop - printf("\n I've put %i particles on the stack \n",nc); + Int_t ipa = iparticle->GetFirstMother()-1; + Int_t iparent = (ipa > -1) ? pParent[ipa] : -1; +// printf("\n SetTrack %d %d %d %d", i, iparent, kf, trackIt[i]); + gAlice->SetTrack(fTrackIt*trackIt[i] , + iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1.); + pParent[i] = nt; + gAlice->KeepTrack(nt); + } // SetTrack loop + } + } else { + nc = GenerateMB(); } // mb ? + if (nc > 0) { jev+=nc; if (jev >= fNpart || fNpart == -1) { @@ -418,93 +420,68 @@ void AliGenPythia::Generate() fXsection=fPythia->GetPARI(1); } -void AliGenPythia::FinishRun() -{ -// Print x-section summary - fPythia->Pystat(1); -} - -Bool_t AliGenPythia::ParentSelected(Int_t ip) -{ -// True if particle is in list of parent particles to be selected - for (Int_t i=0; i<5; i++) - { - if (fParentSelect[i]==ip) return kTRUE; - } - return kFALSE; -} - -Bool_t AliGenPythia::ChildSelected(Int_t ip) +Int_t AliGenPythia::GenerateMB() { -// True if particle is in list of decay products to be selected - if (fForceDecay == kAll) return kTRUE; + Int_t i, kf, nt, iparent; + Int_t nc = 0; + Float_t p[3]; + Float_t polar[3] = {0,0,0}; + Float_t origin[3] = {0,0,0}; + Float_t origin0[3] = {0,0,0}; +// converts from mm/c to s + const Float_t kconv=0.001/2.999792458e8; - for (Int_t i=0; i<5; i++) - { - if (fChildSelect[i]==ip) return kTRUE; - } - return kFALSE; -} + Int_t np = fParticles->GetEntriesFast(); + Int_t* pParent = new Int_t[np]; + for (i=0; i< np-1; i++) pParent[i] = -1; + + for (i = 0; iAt(i); + kf = CheckPDGCode(iparticle->GetPdgCode()); + Int_t ks = iparticle->GetStatusCode(); + Int_t km = iparticle->GetFirstMother(); +// printf("\n Particle: %d %d %d", i, kf, ks); + + if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) || + (ks != 1) || + (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) { + nc++; + if (ks == 1) trackIt = 1; + Int_t ipa = iparticle->GetFirstMother()-1; -Bool_t AliGenPythia::KinematicSelection(TParticle *particle) -{ -// Perform kinematic selection - Float_t px=particle->Px(); - Float_t py=particle->Py(); - Float_t pz=particle->Pz(); - Float_t e=particle->Energy(); + iparent = (ipa > -1) ? pParent[ipa] : -1; // -// transverse momentum cut - Float_t pt=TMath::Sqrt(px*px+py*py); - if (pt > fPtMax || pt < fPtMin) - { -// printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax); - return kFALSE; - } -// -// momentum cut - Float_t p=TMath::Sqrt(px*px+py*py+pz*pz); - if (p > fPMax || p < fPMin) - { -// printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax); - return kFALSE; - } +// store track information + 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.; + Float_t tof=kconv*iparticle->T(); + gAlice->SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar, + tof, kPPrimary, nt); + gAlice->KeepTrack(nt); + pParent[i] = nt; + } // select particle + } // particle loop + + if (pParent) delete[] pParent; -// -// theta cut - Float_t theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz))); - if (theta > fThetaMax || theta < fThetaMin) - { -// printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax); - return kFALSE; - } - -// -// rapidity cut - if ( (e-pz)<=0 || (e+pz)<=0 ) { - return kFALSE; - } - else { - Float_t y = 0.5*TMath::Log((e+pz)/(e-pz)); - if (y > fYMax || y < fYMin) - { -// printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax); - return kFALSE; - } - } + printf("\n I've put %i particles on the stack \n",nc); + MakeHeader(); + return nc; +} -// -// phi cut - Float_t phi=Float_t(TMath::ATan2(Double_t(py),Double_t(px))); - if (phi > fPhiMax || phi < fPhiMin) - { -// printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax); - return kFALSE; - } - return kTRUE; +void AliGenPythia::FinishRun() +{ +// Print x-section summary + fPythia->Pystat(1); } + void AliGenPythia::AdjustWeights() { // Adjust the weights after generation of all events @@ -516,40 +493,6 @@ void AliGenPythia::AdjustWeights() part->SetWeight(part->GetWeight()*fKineBias); } } - -Int_t AliGenPythia::CheckPDGCode(Int_t pdgcode) -{ -// -// If the particle is in a diffractive state, then take action accordingly - switch (pdgcode) { - case 91: - return 92; - case 110: - //rho_diff0 -- difficult to translate, return rho0 - return 113; - case 210: - //pi_diffr+ -- change to pi+ - return 211; - case 220: - //omega_di0 -- change to omega0 - return 223; - case 330: - //phi_diff0 -- return phi0 - return 333; - case 440: - //J/psi_di0 -- return J/psi - return 443; - case 2110: - //n_diffr -- return neutron - return 2112; - case 2210: - //p_diffr+ -- return proton - return 2212; - } - //non diffractive state -- return code unchanged - return pdgcode; -} - void AliGenPythia::SetNuclei(Int_t a1, Int_t a2) { @@ -557,6 +500,15 @@ void AliGenPythia::SetNuclei(Int_t a1, Int_t a2) fNucA1 = a1; fNucA2 = a2; } + + +void AliGenPythia::MakeHeader() +{ +// Builds the event header, to be called after each event + AliGenEventHeader* header = new AliGenPythiaEventHeader("Pythia"); + ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1)); + gAlice->SetGenEventHeader(header); +} AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs) diff --git a/EVGEN/AliGenPythia.h b/EVGEN/AliGenPythia.h index 5057608d59b..daf0e3df416 100644 --- a/EVGEN/AliGenPythia.h +++ b/EVGEN/AliGenPythia.h @@ -6,15 +6,13 @@ /* $Id$ */ -#include "AliGenerator.h" +#include "AliGenMC.h" #include "AliPythia.h" -#include "AliDecayer.h" -#include class AliPythia; class TParticle; -class AliGenPythia : public AliGenerator +class AliGenPythia : public AliGenMC { public: AliGenPythia(); @@ -33,27 +31,26 @@ class AliGenPythia : public AliGenerator {fPtHardMin = ptmin; fPtHardMax = ptmax; } // set centre of mass energy virtual void SetEnergyCMS(Float_t energy = 5500) {fEnergyCMS = energy;} - // force decay type - virtual void SetForceDecay(Decay_t decay = kSemiMuonic) {fForceDecay = decay;} // treat protons as inside nuclei virtual void SetNuclei(Int_t a1, Int_t a2); // get cross section of process virtual Float_t GetXsection() {return fXsection;} - // Check PDG code - virtual Int_t CheckPDGCode(Int_t pdgcode); virtual void FinishRun(); // Assignment Operator AliGenPythia & operator=(const AliGenPythia & rhs); + private: + Int_t GenerateMB(); + virtual void MakeHeader(); protected: + TClonesArray* fParticles; // Particle List + Process_t fProcess; // Process type StrucFunc_t fStrucFunc; // Structure Function - Decay_t fForceDecay; // Decay channel are forced Float_t fEnergyCMS; // Centre of mass energy Float_t fKineBias; // Bias from kinematic selection Int_t fTrials; // Number of trials - TArrayI fParentSelect; // Parent particles to be selected - TArrayI fChildSelect; // Decay products to be selected + Int_t fFlavorSelect; // Heavy Flavor Selection Float_t fXsection; // Cross-section AliPythia *fPythia; //! Pythia Float_t fPtHardMin; // lower pT-hard cut @@ -61,16 +58,10 @@ class AliGenPythia : public AliGenerator Int_t fNucA1; // mass number nucleus side 1 Int_t fNucA2; // mass number nucleus side 2 Bool_t fFullEvent; // Write Full event if true - AliDecayer *fDecayer; // ! pointer to the decayer instance + AliDecayer *fDecayer; //! Pointer to the decayer instance Int_t fDebugEventFirst; // First event to debug Int_t fDebugEventLast; // Last event to debug private: - // check if particle is selected as parent particle - Bool_t ParentSelected(Int_t ip); - // check if particle is selected as child particle - Bool_t ChildSelected(Int_t ip); - // all kinematic selection cuts go here - Bool_t KinematicSelection(TParticle *particle); // adjust the weight from kinematic cuts void AdjustWeights();