1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 Revision 1.6 2000/09/11 13:23:37 morsch
19 Write last seed to file (fortran lun 50) and reed back from same lun using calls to
20 luget_hijing and luset_hijing.
22 Revision 1.5 2000/09/07 16:55:40 morsch
23 fHijing->Initialize(); after change of parameters. (Dmitri Yurevitch Peressounko)
25 Revision 1.4 2000/07/11 18:24:56 fca
26 Coding convention corrections + few minor bug fixes
28 Revision 1.3 2000/06/30 12:08:36 morsch
29 In member data: char* replaced by TString, Init takes care of resizing the strings to
30 8 characters required by Hijing.
32 Revision 1.2 2000/06/15 14:15:05 morsch
33 Add possibility for heavy flavor selection: charm and beauty.
35 Revision 1.1 2000/06/09 20:47:27 morsch
36 AliGenerator interface class to HIJING using THijing (test version)
40 #include "AliGenHijing.h"
41 #include "AliGenHijingEventHeader.h"
46 #include <TParticle.h>
49 ClassImp(AliGenHijing)
51 AliGenHijing::AliGenHijing()
57 AliGenHijing::AliGenHijing(Int_t npart)
60 // Default PbPb collisions at 5. 5 TeV
63 SetImpactParameterRange();
76 AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
82 AliGenHijing::~AliGenHijing()
87 void AliGenHijing::Init()
92 fProjectile.Resize(8);
94 SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget,
95 fAProjectile, fZProjectile, fATarget, fZTarget,
96 fMinImpactParam, fMaxImpactParam));
98 fHijing=(THijing*) fgMCEvGen;
100 fHijing->SetIHPR2(3, fTrigger);
101 fHijing->SetIHPR2(4, fQuench);
102 fHijing->SetIHPR2(6, fShadowing);
103 fHijing->SetIHPR2(12, fDecaysOff);
104 fHijing->SetIHPR2(21, fKeep);
105 fHijing->Rluset(50,0);
106 fHijing->Initialize();
110 if (fEvaluate) EvaluateCrossSections();
113 void AliGenHijing::Generate()
115 // Generate one event
117 Float_t polar[3] = {0,0,0};
118 Float_t origin[3]= {0,0,0};
119 Float_t origin0[3]= {0,0,0};
120 Float_t p[3], random[6];
123 static TClonesArray *particles;
124 // converts from mm/c to s
125 const Float_t kconv=0.001/2.999792458e8;
129 Int_t j, kf, ks, imo;
132 if(!particles) particles=new TClonesArray("TParticle",10000);
135 for (j=0;j<3;j++) origin0[j]=fOrigin[j];
136 if(fVertexSmear==kPerEvent) {
139 origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
140 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
141 // fHijing->SetMSTP(151,0);
143 } else if (fVertexSmear==kPerTrack) {
144 // fHijing->SetMSTP(151,0);
146 // fHijing->SetPARP(151+j, fOsigma[j]*10.);
152 fHijing->GenerateEvent();
154 fHijing->ImportParticles(particles,"All");
155 Int_t np = particles->GetEntriesFast();
156 printf("\n **************************************************%d\n",np);
158 if (np == 0 ) continue;
160 Int_t * newPos = new Int_t[np];
161 for (i = 0; i<np; i++) *(newPos+i)=i;
163 for (i = 0; i<np; i++) {
164 TParticle * iparticle = (TParticle *) particles->At(i);
166 Bool_t hasMother = (iparticle->GetFirstMother() >=0);
167 Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
168 Bool_t selected = kTRUE;
169 Bool_t hasSelectedDaughters = kFALSE;
172 kf = iparticle->GetPdgCode();
173 if (!fSelectAll) selected = KinematicSelection(iparticle)&&SelectFlavor(kf);
174 if (hasDaughter && !selected) hasSelectedDaughters = DaughtersSelection(iparticle, particles);
176 // Put particle on the stack if it is either selected or it is the mother of at least one seleted particle
179 if (selected || hasSelectedDaughters) {
181 ks = iparticle->GetStatusCode();
182 p[0]=iparticle->Px();
183 p[1]=iparticle->Py();
184 p[2]=iparticle->Pz();
185 origin[0]=origin0[0]+iparticle->Vx()/10;
186 origin[1]=origin0[1]+iparticle->Vy()/10;
187 origin[2]=origin0[2]+iparticle->Vz()/10;
188 tof=kconv*iparticle->T();
191 imo=iparticle->GetFirstMother();
195 // printf("\n selected iparent %d %d %d \n",i, kf, imo);
197 gAlice->SetTrack(0,imo,kf,p,origin,polar,
200 gAlice->SetTrack(fTrackIt,imo,kf,p,origin,polar,
207 printf("\n I've put %i particles on the stack \n",nc);
210 if (jev >= fNpart || fNpart == -1) {
211 fKineBias=Float_t(fNpart)/Float_t(fTrials);
212 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
217 fHijing->Rluget(50,-1);
220 Bool_t AliGenHijing::KinematicSelection(TParticle *particle)
222 // Perform kinematic selection
223 Float_t px=particle->Px();
224 Float_t py=particle->Py();
225 Float_t pz=particle->Pz();
226 Float_t e=particle->Energy();
229 // transverse momentum cut
230 Float_t pt=TMath::Sqrt(px*px+py*py);
231 if (pt > fPtMax || pt < fPtMin)
233 // printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
238 Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
239 if (p > fPMax || p < fPMin)
241 // printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
247 Float_t theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
248 if (theta > fThetaMax || theta < fThetaMin)
251 // printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
257 Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
258 if (y > fYMax || y < fYMin)
260 // printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
266 Float_t phi=Float_t(TMath::ATan2(Double_t(py),Double_t(px)));
267 if (phi > fPhiMax || phi < fPhiMin)
269 // printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
276 void AliGenHijing::KeepFullEvent()
281 void AliGenHijing::EvaluateCrossSections()
283 // Glauber Calculation of geometrical x-section
285 Float_t xTot=0.; // barn
286 Float_t xTotHard=0.; // barn
287 Float_t xPart=0.; // barn
288 Float_t xPartHard=0.; // barn
289 Float_t sigmaHard=0.1; // mbarn
291 Float_t bMax=fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
292 const Float_t kdib=0.2;
293 Int_t kMax=Int_t((bMax-bMin)/kdib)+1;
296 printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
297 printf("\n Target Radius (fm): %f \n",fHijing->GetHIPR1(35));
301 for (i=0; i<kMax; i++)
303 Float_t xb=bMin+i*kdib;
305 ov=fHijing->Profile(xb);
306 Float_t gb = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
307 Float_t gbh = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
310 if (xb > fMinImpactParam && xb < fMaxImpactParam)
316 if ((xTot-oldvalue)/oldvalue<0.0001) break;
318 printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
319 printf("\n Hard cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
321 printf("\n Total cross section (barn): %f \n",xTot);
322 printf("\n Hard cross section (barn): %f \n \n",xTotHard);
323 printf("\n Partial cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
324 printf("\n Partial hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
327 Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
330 // Looks recursively if one of the daughters has been selected
332 // printf("\n Consider daughters %d:",iparticle->GetPdgCode());
336 Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
337 Bool_t selected=kFALSE;
339 imin=iparticle->GetFirstDaughter();
340 imax=iparticle->GetLastDaughter();
341 for (i=imin; i<= imax; i++){
342 TParticle * jparticle = (TParticle *) particles->At(i);
343 Int_t ip=jparticle->GetPdgCode();
344 if (KinematicSelection(jparticle)&&SelectFlavor(ip)) {selected=kTRUE; break;}
345 if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
355 Bool_t AliGenHijing::SelectFlavor(Int_t pid)
357 // Select flavor of particle
359 // 4: charm and beauty
361 if (fFlavor == 0) return kTRUE;
363 Int_t ifl=TMath::Abs(pid/100);
364 if (ifl > 10) ifl/=10;
365 return ((fFlavor==4 && (ifl==4 || ifl==5)) ||
366 (fFlavor==5 && ifl==5));
370 void AliGenHijing::MakeHeader()
372 // Builds the event header, to be called after each event
373 AliGenHijingEventHeader* header = new AliGenHijingEventHeader("Hijing");
374 // header->SetDate(date);
375 // header->SetRunNumber(run);
376 // header->SetEventNumber(event);
377 header->SetNProduced(fHijing->GetNATT());
378 header->SetImpactParameter(fHijing->GetHINT1(19));
379 header->SetTotalEnergy(fHijing->GetEATT());
380 header->SetHardScatters(fHijing->GetJATT());
381 header->SetParticipants(fHijing->GetNP(), fHijing->GetNT());
382 header->SetCollisions(fHijing->GetN0(),
388 AliGenHijing& AliGenHijing::operator=(const AliGenHijing& rhs)
390 // Assignment operator