X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EVGEN%2FAliGenHijing.cxx;h=17b37c277e683d25ace4c7980c197c26dd6f4076;hb=ff0f2bbda7b5113976cee1be92e458f796bc795f;hp=990764c0a657524e05a849a153b38fa5ca590bc7;hpb=7498b3ad12680aeae30cc037737c9f0e7de876d1;p=u%2Fmrichter%2FAliRoot.git diff --git a/EVGEN/AliGenHijing.cxx b/EVGEN/AliGenHijing.cxx index 990764c0a65..17b37c277e6 100644 --- a/EVGEN/AliGenHijing.cxx +++ b/EVGEN/AliGenHijing.cxx @@ -15,6 +15,82 @@ /* $Log$ +Revision 1.40 2002/02/12 11:05:53 morsch +Get daughter indices right. + +Revision 1.39 2002/02/12 09:16:39 morsch +Correction in SelectFlavor() + +Revision 1.38 2002/02/12 08:53:21 morsch +SetNoGammas can be used to inhibit writing of gammas and pi0. + +Revision 1.37 2002/02/08 16:50:50 morsch +Add name and title in constructor. + +Revision 1.36 2002/01/31 20:17:55 morsch +Allow for triggered jets with simplified topology: Exact pT, back-to-back + +Revision 1.35 2001/12/13 07:56:25 hristov +Set pointers to zero in the default constructor + +Revision 1.34 2001/12/11 16:55:42 morsch +Correct initialization for jet phi-range. + +Revision 1.33 2001/12/05 10:18:51 morsch +Possibility of kinematic biasing of jet phi range. (J. Klay) + +Revision 1.32 2001/11/28 13:51:11 morsch +Introduce kinematic biasing (etamin, etamax) of jet trigger. Bookkeeping +(number of trials) done in AliGenHijingEventHeader. + +Revision 1.31 2001/11/06 12:30:34 morsch +Add Boost() method to boost all particles to LHC lab frame. Needed for asymmetric collision systems. + +Revision 1.30 2001/10/21 18:35:56 hristov +Several pointers were set to zero in the default constructors to avoid memory management problems + +Revision 1.29 2001/10/15 08:12:24 morsch +- Vertex smearing with truncated gaussian. +- Store triggered jet info before and after final state radiation into mc-heade + +Revision 1.28 2001/10/08 11:55:25 morsch +Store 4-momenta of trigegred jets in event header. +Possibility to switch of initial and final state radiation. + +Revision 1.27 2001/10/08 07:13:14 morsch +Add setter for minimum transverse momentum of triggered jet. + +Revision 1.26 2001/10/04 08:12:24 morsch +Redefinition of stable condition. + +Revision 1.25 2001/07/27 17:09:36 morsch +Use local SetTrack, KeepTrack and SetHighWaterMark methods +to delegate either to local stack or to stack owned by AliRun. +(Piotr Skowronski, A.M.) + +Revision 1.24 2001/07/20 09:34:56 morsch +Count the number of spectator neutrons and protons and add information +to the event header. (Chiara Oppedisano) + +Revision 1.23 2001/07/13 17:30:22 morsch +Derive from AliGenMC. + +Revision 1.22 2001/06/11 13:09:23 morsch +- Store cross-Section and number of binary collisions as a function of impact parameter +- Pass AliGenHijingEventHeader to gAlice. + +Revision 1.21 2001/02/28 17:35:24 morsch +Consider elastic interactions (ks = 1 and ks = 11) as spectator (Chiara Oppedisano) + +Revision 1.20 2001/02/14 15:50:40 hristov +The last particle in event marked using SetHighWaterMark + +Revision 1.19 2000/12/21 16:24:06 morsch +Coding convention clean-up + +Revision 1.18 2000/12/06 17:46:30 morsch +Avoid random numbers 1 and 0. + Revision 1.17 2000/12/04 11:22:03 morsch Init of sRandom as in 1.15 @@ -73,44 +149,76 @@ AliGenerator interface class to HIJING using THijing (test version) */ + + +// Generator using HIJING as an external generator +// The main HIJING options are accessable for the user through this interface. +// Uses the THijing implementation of TGenerator. +// +// andreas.morsch@cern.ch + #include "AliGenHijing.h" #include "AliGenHijingEventHeader.h" #include "AliRun.h" +#include "AliPDG.h" #include #include #include +#include +#include ClassImp(AliGenHijing) AliGenHijing::AliGenHijing() - :AliGenerator() + :AliGenMC() { // Constructor + fHijing = 0; + fDsigmaDb = 0; + fDnDb = 0; } AliGenHijing::AliGenHijing(Int_t npart) - :AliGenerator(npart) + :AliGenMC(npart) { // Default PbPb collisions at 5. 5 TeV // + fName = "Hijing"; + fTitle= "Particle Generator using HIJING"; + SetEnergyCMS(); SetImpactParameterRange(); SetTarget(); SetProjectile(); - fKeep=0; - fQuench=1; - fShadowing=1; - fTrigger=0; - fDecaysOff=1; - fEvaluate=0; - fSelectAll=0; - fFlavor=0; - fSpectators=1; + SetBoostLHC(); + SetJetEtaRange(); + SetJetPhiRange(); + + fKeep = 0; + fQuench = 1; + fShadowing = 1; + fTrigger = 0; + fDecaysOff = 1; + fEvaluate = 0; + fSelectAll = 0; + fFlavor = 0; + fSpectators = 1; + fDsigmaDb = 0; + fDnDb = 0; + fPtMinJet = -2.5; + fRadiation = 1; + fEventVertex.Set(3); +// + SetSimpleJets(); + SetNoGammas(); + // // Set random number generator sRandom = fRandom; + fHijing = 0; + } AliGenHijing::AliGenHijing(const AliGenHijing & Hijing) @@ -122,6 +230,8 @@ AliGenHijing::AliGenHijing(const AliGenHijing & Hijing) AliGenHijing::~AliGenHijing() { // Destructor + if ( fDsigmaDb) delete fDsigmaDb; + if ( fDnDb) delete fDnDb; } void AliGenHijing::Init() @@ -136,16 +246,51 @@ void AliGenHijing::Init() fMinImpactParam, fMaxImpactParam)); fHijing=(THijing*) fgMCEvGen; - + fHijing->SetIHPR2(2, fRadiation); fHijing->SetIHPR2(3, fTrigger); - fHijing->SetIHPR2(4, fQuench); fHijing->SetIHPR2(6, fShadowing); fHijing->SetIHPR2(12, fDecaysOff); fHijing->SetIHPR2(21, fKeep); - fHijing->Rluset(50,0); - fHijing->Initialize(); - + fHijing->SetHIPR1(10, fPtMinJet); + fHijing->SetHIPR1(50, fSimpleJet); +// +// Quenching +// +// +// fQuench = 0: no quenching +// fQuench = 1: hijing default +// fQuench = 2: new LHC parameters for HIPR1(11) and HIPR1(14) +// fQuench = 3: new RHIC parameters for HIPR1(11) and HIPR1(14) +// fQuench = 4: new LHC parameters with log(e) dependence +// fQuench = 5: new RHIC parameters with log(e) dependence + fHijing->SetIHPR2(50, 0); + if (fQuench > 0) + fHijing->SetIHPR2(4, 1); + else + fHijing->SetIHPR2(4, 0); +// New LHC parameters from Xin-Nian Wang + if (fQuench == 2) { + fHijing->SetHIPR1(14, 1.1); + fHijing->SetHIPR1(11, 3.7); + } else if (fQuench == 3) { + fHijing->SetHIPR1(14, 0.20); + fHijing->SetHIPR1(11, 2.5); + } else if (fQuench == 4) { + fHijing->SetIHPR2(50, 1); + fHijing->SetHIPR1(14, 4.*0.34); + fHijing->SetHIPR1(11, 3.7); + } else if (fQuench == 5) { + fHijing->SetIHPR2(50, 1); + fHijing->SetHIPR1(14, 0.34); + fHijing->SetHIPR1(11, 2.5); + } + + +// +// Initialize Hijing +// + fHijing->Initialize(); // if (fEvaluate) EvaluateCrossSections(); // @@ -157,220 +302,187 @@ void AliGenHijing::Generate() { // Generate one event - 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 tof; + 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 tof; - static TClonesArray *particles; + static TClonesArray *particles; // converts from mm/c to s - const Float_t kconv=0.001/2.999792458e8; + const Float_t kconv = 0.001/2.999792458e8; // - Int_t nt=0; - Int_t jev=0; - Int_t j, kf, ks, imo; - kf=0; + Int_t nt = 0; + Int_t jev = 0; + Int_t j, kf, ks, imo; + kf = 0; - if(!particles) particles=new TClonesArray("TParticle",10000); + if(!particles) particles = new TClonesArray("TParticle",10000); - fTrials=0; - 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])); + 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]; + } else if (fVertexSmear == kPerTrack) { // fHijing->SetMSTP(151,0); - } - } else if (fVertexSmear==kPerTrack) { -// fHijing->SetMSTP(151,0); - for (j=0;j<3;j++) { -// fHijing->SetPARP(151+j, fOsigma[j]*10.); - } - } - while(1) - { - - fHijing->GenerateEvent(); - fTrials++; - fHijing->ImportParticles(particles,"All"); - Int_t np = particles->GetEntriesFast(); - printf("\n **************************************************%d\n",np); - Int_t nc=0; - if (np == 0 ) continue; - Int_t i; - Int_t * newPos = new Int_t[np]; - - for (i = 0; iSetPARP(151+j, fOsigma[j]*10.); + } + } + while(1) + { +// Generate one event +// -------------------------------------------------------------------------- + fSpecn = 0; + fSpecp = 0; +// -------------------------------------------------------------------------- + fHijing->GenerateEvent(); + fTrials++; + if (fTrigger != kNoTrigger) { + if (!CheckTrigger()) continue; + } + + fHijing->ImportParticles(particles,"All"); + if (fLHC) Boost(particles); + + + Int_t np = particles->GetEntriesFast(); + printf("\n **************************************************%d\n",np); + Int_t nc = 0; + if (np == 0 ) continue; + Int_t i; + Int_t* newPos = new Int_t[np]; + Int_t* pSelected = new Int_t[np]; + + for (i = 0; i < np; i++) { + newPos[i] = i; + pSelected[i] = 0; + } + +// Get event vertex // -// First write parent particles + TParticle * iparticle = (TParticle *) particles->At(0); + fEventVertex[0] = origin0[0]; + fEventVertex[1] = origin0[1]; + fEventVertex[2] = origin0[2]; + +// +// First select parent particles // - for (i = 0; iAt(i); + for (i = 0; i < np; i++) { + iparticle = (TParticle *) particles->At(i); // Is this a parent particle ? - if (Stable(iparticle)) continue; + if (Stable(iparticle)) continue; // - Bool_t hasMother = (iparticle->GetFirstMother() >=0); - Bool_t selected = kTRUE; - Bool_t hasSelectedDaughters = kFALSE; - + Bool_t selected = kTRUE; + Bool_t hasSelectedDaughters = kFALSE; + + + kf = iparticle->GetPdgCode(); + ks = iparticle->GetStatusCode(); + if (kf == 92) continue; - kf = iparticle->GetPdgCode(); - ks = iparticle->GetStatusCode(); - if (kf == 92) continue; - - if (!fSelectAll) selected = KinematicSelection(iparticle)&&SelectFlavor(kf); - hasSelectedDaughters = DaughtersSelection(iparticle, particles); -// -// Put particle on the stack if it is either selected or it is the mother of at least one seleted particle -// - if (selected || hasSelectedDaughters) { - nc++; - 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; - tof=kconv*iparticle->T(); - imo=-1; - if (hasMother) { - imo=iparticle->GetFirstMother(); - TParticle* mother= (TParticle *) particles->At(imo); - imo = (mother->GetPdgCode() != 92) ? imo=*(newPos+imo) : -1; - } -// Put particle on the stack ... -// printf("\n set track mother: %d %d %d %d %d %d ",i,imo, kf, nt+1, selected, hasSelectedDaughters); - - gAlice->SetTrack(0,imo,kf,p,origin,polar, - tof,kPPrimary,nt); -// ... and keep it there - gAlice->KeepTrack(nt); -// - *(newPos+i)=nt; - } // selected - } // particle loop parents -// -// Now write the final state particles -// - - for (i = 0; iAt(i); -// Is this a final state particle ? - if (!Stable(iparticle)) continue; -// - Bool_t hasMother = (iparticle->GetFirstMother() >=0); - Bool_t selected = kTRUE; - kf = iparticle->GetPdgCode(); - ks = iparticle->GetStatusCode(); - if (!fSelectAll) { - selected = KinematicSelection(iparticle)&&SelectFlavor(kf); - if (!fSpectators && selected) selected = (ks != 0 && ks != 10); - } -// -// Put particle on the stack if selected + if (!fSelectAll) selected = KinematicSelection(iparticle, 0) && + SelectFlavor(kf); + hasSelectedDaughters = DaughtersSelection(iparticle, particles); // - if (selected) { - nc++; - 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; - tof=kconv*iparticle->T(); - imo=-1; - - if (hasMother) { - imo=iparticle->GetFirstMother(); - TParticle* mother= (TParticle *) particles->At(imo); - imo = (mother->GetPdgCode() != 92) ? imo=*(newPos+imo) : -1; - } -// Put particle on the stack - gAlice->SetTrack(fTrackIt,imo,kf,p,origin,polar, - tof,kPNoProcess,nt); -// tof,"Secondary",nt); - -// printf("\n set track final: %d %d %d",imo, kf, nt); - gAlice->KeepTrack(nt); - *(newPos+i)=nt; - } // selected - } // particle loop final state - - delete newPos; - - printf("\n I've put %i particles on the stack \n",nc); - if (nc > 0) { - jev+=nc; - if (jev >= fNpart || fNpart == -1) { - fKineBias=Float_t(fNpart)/Float_t(fTrials); - printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev); - break; - } - } - } // event loop - fHijing->Rluget(50,-1); -} - -Bool_t AliGenHijing::KinematicSelection(TParticle *particle) -{ -// Perform kinematic selection - Double_t px=particle->Px(); - Double_t py=particle->Py(); - Double_t pz=particle->Pz(); - Double_t e=particle->Energy(); - +// Put particle on the stack if it is either selected or +// it is the mother of at least one seleted particle // -// transverse momentum cut - Double_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; - } + if (selected || hasSelectedDaughters) { + nc++; + pSelected[i] = 1; + } // selected + } // particle loop parents // -// momentum cut - Double_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; - } - +// Now select the final state particles // -// theta cut - Double_t theta = Double_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; - } + for (i = 0; iAt(i); +// Is this a final state particle ? + if (!Stable(iparticle)) continue; + + Bool_t selected = kTRUE; + kf = iparticle->GetPdgCode(); + ks = iparticle->GetStatusCode(); +// -------------------------------------------------------------------------- +// Count spectator neutrons and protons + if(ks == 0 || ks == 1 || ks == 10 || ks == 11){ + if(kf == kNeutron) fSpecn += 1; + if(kf == kProton) fSpecp += 1; + } +// -------------------------------------------------------------------------- +// + if (!fSelectAll) { + selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf); + if (!fSpectators && selected) selected = (ks != 0 && ks != 1 && ks != 10 + && ks != 11); + } // -// rapidity cut - Double_t y; - if(e<=pz) y = 99; - else if (e<=-pz) y = -99; - else 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; - } - +// Put particle on the stack if selected // -// phi cut - Double_t phi=Double_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; + if (selected) { + nc++; + pSelected[i] = 1; + } // selected + } // particle loop final state +// +// Write particles to stack +// + for (i = 0; iAt(i); + Bool_t hasMother = (iparticle->GetFirstMother() >=0); + Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0); + + if (pSelected[i]) { + kf = iparticle->GetPdgCode(); + 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; + tof = kconv*iparticle->T(); + imo = -1; + TParticle* mother = 0; + if (hasMother) { + imo = iparticle->GetFirstMother(); + mother = (TParticle *) particles->At(imo); + imo = (mother->GetPdgCode() != 92) ? imo = newPos[imo] : -1; + } // if has mother + Bool_t tFlag = (fTrackIt && !hasDaughter); + SetTrack(tFlag,imo,kf,p,origin,polar, + tof,kPNoProcess,nt); + 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) { + fKineBias = Float_t(fNpart)/Float_t(fTrials); + printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev); + break; + } + } + } // event loop + MakeHeader(); + SetHighWaterMark(nt); } void AliGenHijing::KeepFullEvent() @@ -382,46 +494,64 @@ void AliGenHijing::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=fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35); - const Float_t kdib=0.2; - Int_t kMax=Int_t((bMax-bMin)/kdib)+1; + 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 = fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35); + const Float_t kdib = 0.2; + Int_t kMax = Int_t((bMax-bMin)/kdib)+1; printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34)); printf("\n Target Radius (fm): %f \n",fHijing->GetHIPR1(35)); Int_t i; - Float_t oldvalue=0.; + 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; iProfile(xb); - Float_t gb = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov)); - Float_t gbh = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov; + Float_t gb = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov)); + Float_t gbh = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov; xTot+=gb; - xTotHard+=gbh; + xTotHard += gbh; if (xb > fMinImpactParam && xb < fMaxImpactParam) { - xPart+=gb; - xPartHard+=gbh; + xPart += gb; + xPartHard += gbh; } if(oldvalue) if ((xTot-oldvalue)/oldvalue<0.0001) break; - oldvalue=xTot; + 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 AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* particles) @@ -430,18 +560,18 @@ Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* part // Looks recursively if one of the daughters has been selected // // printf("\n Consider daughters %d:",iparticle->GetPdgCode()); - Int_t imin=-1; - Int_t imax=-1; + Int_t imin = -1; + Int_t imax = -1; Int_t i; - Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0); - Bool_t selected=kFALSE; + Bool_t hasDaughters = (iparticle->GetFirstDaughter() >=0); + Bool_t selected = kFALSE; if (hasDaughters) { - imin=iparticle->GetFirstDaughter(); - imax=iparticle->GetLastDaughter(); - for (i=imin; i<= imax; i++){ - TParticle * jparticle = (TParticle *) particles->At(i); - Int_t ip=jparticle->GetPdgCode(); - if (KinematicSelection(jparticle)&&SelectFlavor(ip)) { + imin = iparticle->GetFirstDaughter(); + imax = iparticle->GetLastDaughter(); + for (i = imin; i <= imax; i++){ + TParticle * jparticle = (TParticle *) particles->At(i); + Int_t ip = jparticle->GetPdgCode(); + if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) { selected=kTRUE; break; } if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; } @@ -449,7 +579,6 @@ Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* part } else { return kFALSE; } - return selected; } @@ -460,19 +589,29 @@ Bool_t AliGenHijing::SelectFlavor(Int_t pid) // 0: all // 4: charm and beauty // 5: beauty - if (fFlavor == 0) return kTRUE; + Bool_t res = 0; - Int_t ifl=TMath::Abs(pid/100); - if (ifl > 10) ifl/=10; - return (fFlavor == ifl); + if (fFlavor == 0) { + res = kTRUE; + } else { + Int_t ifl = TMath::Abs(pid/100); + if (ifl > 10) ifl/=10; + res = (fFlavor == ifl); + } +// +// This part if gamma writing is inhibited + if (fNoGammas) + res = res && (pid != kGamma && pid != kPi0); +// + return res; } Bool_t AliGenHijing::Stable(TParticle* particle) { - Int_t kf = TMath::Abs(particle->GetPdgCode()); +// Return true for a stable particle +// - if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122)) - + if (particle->GetFirstDaughter() < 0 ) { return kTRUE; } else { @@ -480,24 +619,123 @@ Bool_t AliGenHijing::Stable(TParticle* particle) } } + +void AliGenHijing::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*) particles->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 AliGenHijing::MakeHeader() { // Builds the event header, to be called after each event - AliGenHijingEventHeader* header = new AliGenHijingEventHeader("Hijing"); -// header->SetDate(date); -// header->SetRunNumber(run); -// header->SetEventNumber(event); - header->SetNProduced(fHijing->GetNATT()); - header->SetImpactParameter(fHijing->GetHINT1(19)); - header->SetTotalEnergy(fHijing->GetEATT()); - header->SetHardScatters(fHijing->GetJATT()); - header->SetParticipants(fHijing->GetNP(), fHijing->GetNT()); - header->SetCollisions(fHijing->GetN0(), - fHijing->GetN01(), - fHijing->GetN10(), - fHijing->GetN11()); + AliGenEventHeader* header = new AliGenHijingEventHeader("Hijing"); + ((AliGenHijingEventHeader*) header)->SetNProduced(fHijing->GetNATT()); + ((AliGenHijingEventHeader*) header)->SetImpactParameter(fHijing->GetHINT1(19)); + ((AliGenHijingEventHeader*) header)->SetTotalEnergy(fHijing->GetEATT()); + ((AliGenHijingEventHeader*) header)->SetHardScatters(fHijing->GetJATT()); + ((AliGenHijingEventHeader*) header)->SetParticipants(fHijing->GetNP(), fHijing->GetNT()); + ((AliGenHijingEventHeader*) header)->SetCollisions(fHijing->GetN0(), + fHijing->GetN01(), + fHijing->GetN10(), + fHijing->GetN11()); + ((AliGenHijingEventHeader*) header)->SetSpectators(fSpecn, fSpecp); + +// 4-momentum vectors of the triggered jets. +// +// Before final state gluon radiation. + TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21), + fHijing->GetHINT1(22), + fHijing->GetHINT1(23), + fHijing->GetHINT1(24)); + + TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31), + fHijing->GetHINT1(32), + fHijing->GetHINT1(33), + fHijing->GetHINT1(34)); +// After final state gluon radiation. + TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26), + fHijing->GetHINT1(27), + fHijing->GetHINT1(28), + fHijing->GetHINT1(29)); + + TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36), + fHijing->GetHINT1(37), + fHijing->GetHINT1(38), + fHijing->GetHINT1(39)); + ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4); +// Bookkeeping for kinematic bias + ((AliGenHijingEventHeader*) header)->SetTrials(fTrials); +// Event Vertex + header->SetPrimaryVertex(fEventVertex); + gAlice->SetGenEventHeader(header); } +Bool_t AliGenHijing::CheckTrigger() +{ +// Check the kinematic trigger condition +// + TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(26), + fHijing->GetHINT1(27), + fHijing->GetHINT1(28), + fHijing->GetHINT1(29)); + + TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(36), + fHijing->GetHINT1(37), + fHijing->GetHINT1(38), + fHijing->GetHINT1(39)); + Double_t eta1 = jet1->Eta(); + Double_t eta2 = jet2->Eta(); + Double_t phi1 = jet1->Phi(); + Double_t phi2 = jet2->Phi(); + Bool_t triggered = kFALSE; +// printf("\n Trigger: %f %f %f %f", +// fEtaMinJet, fEtaMaxJet, fPhiMinJet, fPhiMaxJet); +// printf("\n Jet1: %f %f", phi1, eta1); +// printf("\n Jet2: %f %f", phi2, eta2); + + + if ( + (eta1 < fEtaMaxJet && eta1 > fEtaMinJet && + phi1 < fPhiMaxJet && phi1 > fPhiMinJet) + || + (eta2 < fEtaMaxJet && eta2 > fEtaMinJet && + phi2 < fPhiMaxJet && phi2 > fPhiMinJet) + ) + triggered = kTRUE; + + return triggered; +} + + + + AliGenHijing& AliGenHijing::operator=(const AliGenHijing& rhs) { // Assignment operator