JYMMY includes.
[u/mrichter/AliRoot.git] / THerwig / AliGenHerwig.cxx
CommitLineData
618e1dc0 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
7cdba479 16/* $Id$ */
618e1dc0 17
18
19// Generator using Herwig as an external generator
20// The main Herwig options are accessable for the user through this interface.
21// Uses the THerwig implementation of TGenerator.
22
23#include "AliGenHerwig.h"
93831bcf 24#include "AliHerwigRndm.h"
618e1dc0 25#include "AliRun.h"
26
27#include <TParticle.h>
28#include "THerwig6.h"
29
30#include "Riostream.h"
5d12ce38 31#include "AliMC.h"
618e1dc0 32
e2054d85 33#include "driver.h"
34
014a9521 35ClassImp(AliGenHerwig)
618e1dc0 36
7cdba479 37
a1d47d99 38 AliGenHerwig::AliGenHerwig() :
39 AliGenMC(),
40 fAutPDF("GRV"),
41 fModPDF(5),
2ad1e1c2 42 fStrucFunc(kCTEQ5L),
a1d47d99 43 fKeep(0),
44 fDecaysOff(1),
45 fTrigger(0),
46 fSelectAll(0),
47 fFlavor(0),
48 fEnergyCMS(14000),
49 fMomentum1(7000),
50 fMomentum2(7000),
51 fKineBias(1),
52 fTrials(0),
53 fXsection(0),
54 fHerwig(0x0),
55 fProcess(0),
2d654757 56 fPtHardMin(0.),
57 fPtRMS(0.),
a1d47d99 58 fMaxPr(10),
59 fMaxErrors(1000),
e2054d85 60 fEnSoft(1),
61 fEv1Pr(0),
2d654757 62 fEv2Pr(0),
63 fFileName(0)
618e1dc0 64{
65// Constructor
66}
67
68AliGenHerwig::AliGenHerwig(Int_t npart)
2d654757 69 :AliGenMC(npart),
70 fAutPDF("GRV"),
71 fModPDF(5),
72 fStrucFunc(kCTEQ5L),
73 fKeep(0),
74 fDecaysOff(1),
75 fTrigger(0),
76 fSelectAll(0),
77 fFlavor(0),
78 fEnergyCMS(14000),
79 fMomentum1(7000),
80 fMomentum2(7000),
81 fKineBias(1),
82 fTrials(0),
83 fXsection(0),
84 fHerwig(0x0),
85 fProcess(0),
86 fPtHardMin(0.),
87 fPtRMS(0.),
88 fMaxPr(10),
89 fMaxErrors(1000),
90 fEnSoft(1),
91 fEv1Pr(0),
92 fEv2Pr(0),
93 fFileName(0)
618e1dc0 94{
618e1dc0 95 SetTarget();
96 SetProjectile();
93831bcf 97 // Set random number generator
98 AliHerwigRndm::SetHerwigRandom(GetRandom());
618e1dc0 99}
100
101AliGenHerwig::AliGenHerwig(const AliGenHerwig & Herwig)
2d654757 102 :AliGenMC(Herwig),
103 fAutPDF("GRV"),
104 fModPDF(5),
105 fStrucFunc(kCTEQ5L),
106 fKeep(0),
107 fDecaysOff(1),
108 fTrigger(0),
109 fSelectAll(0),
110 fFlavor(0),
111 fEnergyCMS(14000),
112 fMomentum1(7000),
113 fMomentum2(7000),
114 fKineBias(1),
115 fTrials(0),
116 fXsection(0),
117 fHerwig(0x0),
118 fProcess(0),
119 fPtHardMin(0.),
120 fPtRMS(0.),
121 fMaxPr(10),
122 fMaxErrors(1000),
123 fEnSoft(1),
124 fEv1Pr(0),
125 fEv2Pr(0),
126 fFileName(0)
618e1dc0 127{
014a9521 128// Copy constructor
129 Herwig.Copy(*this);
618e1dc0 130}
131
132
133AliGenHerwig::~AliGenHerwig()
134{
135// Destructor
136}
137
e2054d85 138void AliGenHerwig::SetEventListRange(Int_t eventFirst, Int_t eventLast)
139{
140 fEv1Pr = ++eventFirst;
141 fEv2Pr = ++eventLast;
142 if ( fEv2Pr == -1 ) fEv2Pr = fEv2Pr;
143}
144
618e1dc0 145void AliGenHerwig::Init()
146{
147// Initialisation
148 fTarget.Resize(8);
149 fProjectile.Resize(8);
150 SetMC(new THerwig6());
d274cd60 151 fHerwig=(THerwig6*) fMCEvGen;
618e1dc0 152 // initialize common blocks
153 fHerwig->Initialize(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
154 // reset parameters according to user needs
155 InitPDF();
156 fHerwig->SetPTMIN(fPtHardMin);
157 fHerwig->SetPTRMS(fPtRMS);
158 fHerwig->SetMAXPR(fMaxPr);
159 fHerwig->SetMAXER(fMaxErrors);
160 fHerwig->SetENSOF(fEnSoft);
e2054d85 161
162 fHerwig->SetEV1PR(fEv1Pr);
163 fHerwig->SetEV2PR(fEv2Pr);
164
165// C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
166// RMASS(1)=0.32
167// RMASS(2)=0.32
168// RMASS(3)=0.5
169// RMASS(4)=1.55
170// RMASS(5)=4.95
171// RMASS(6)=174.3
172// RMASS(13)=0.75
173
174 fHerwig->SetRMASS(4,1.2);
175 fHerwig->SetRMASS(5,4.75);
176
177 if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
178
179 fHerwig->Hwusta("PI0 ");
180
618e1dc0 181 // compute parameter dependent constants
182 fHerwig->PrepareRun();
183}
184
185void AliGenHerwig::InitPDF()
186{
187 switch(fStrucFunc)
188 {
2d654757 189// case kGRVLO:
190// fModPDF=5;
191// fAutPDF="GRV";
192// break;
193// case kGRVHO:
194// fModPDF=6;
195// fAutPDF="GRV";
196// break;
618e1dc0 197 case kGRVLO98:
198 fModPDF=12;
199 fAutPDF="GRV";
200 break;
2d654757 201// case kMRSDminus:
202// fModPDF=31;
203// fAutPDF="MRS";
204// break;
205// case kMRSD0:
206// fModPDF=30;
207// fAutPDF="MRS";
208// break;
209// case kMRSG:
210// fModPDF=41;
211// fAutPDF="MRS";
212// break;
213// case kMRSTcgLO:
214// fModPDF=72;
215// fAutPDF="MRS";
216// break;
618e1dc0 217 case kCTEQ4M:
218 fModPDF=34;
219 fAutPDF="CTEQ";
220 break;
221 case kCTEQ5L:
222 fModPDF=46;
223 fAutPDF="CTEQ";
224 break;
2d654757 225// case kCTEQ5M:
226// fModPDF=48;
227// fAutPDF="CTEQ";
228// break;
618e1dc0 229 default:
230 cerr << "This structure function is not inplemented " << fStrucFunc << endl;
231 break;
232 }
233 fAutPDF.Resize(20);
234 fHerwig->SetMODPDF(1,fModPDF);
235 fHerwig->SetMODPDF(2,fModPDF);
236 fHerwig->SetAUTPDF(1,fAutPDF);
237 fHerwig->SetAUTPDF(2,fAutPDF);
238}
239
240void AliGenHerwig::Generate()
241{
242 // Generate one event
243
244 Float_t polar[3] = {0,0,0};
245 Float_t origin[3]= {0,0,0};
246 Float_t origin0[3]= {0,0,0};
247 Float_t p[4], random[6];
248
249 static TClonesArray *particles;
250 // converts from mm/c to s
251 const Float_t kconv=0.001/2.999792458e8;
252 //
253 Int_t nt=0;
254 Int_t jev=0;
255 Int_t j, kf, ks, imo;
256 kf=0;
257
258 if(!particles) particles=new TClonesArray("TParticle",10000);
259
260 fTrials=0;
261 for (j=0;j<3;j++) origin0[j]=fOrigin[j];
262 if(fVertexSmear==kPerEvent) {
263 Rndm(random,6);
264 for (j=0;j<3;j++) {
265 origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
266 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
267 }
268 }
269
270 while(1)
271 {
272 fHerwig->GenerateEvent();
273 fTrials++;
274 fHerwig->ImportParticles(particles,"All");
275 Int_t np = particles->GetEntriesFast()-1;
276 if (np == 0 ) continue;
277
278 Int_t nc=0;
279
280 Int_t * newPos = new Int_t[np];
281 for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
282
283 for (Int_t i = 0; i<np; i++) {
284 TParticle * iparticle = (TParticle *) particles->At(i);
285 imo = iparticle->GetFirstMother();
286 kf = iparticle->GetPdgCode();
287 ks = iparticle->GetStatusCode();
288 if (ks != 3 &&
289 KinematicSelection(iparticle,0))
290 {
291 nc++;
292 p[0]=iparticle->Px();
293 p[1]=iparticle->Py();
294 p[2]=iparticle->Pz();
295 p[3]=iparticle->Energy();
296 origin[0]=origin0[0]+iparticle->Vx()/10;
297 origin[1]=origin0[1]+iparticle->Vy()/10;
298 origin[2]=origin0[2]+iparticle->Vz()/10;
299 Float_t tof = kconv*iparticle->T();
300 Int_t iparent = (imo > -1) ? newPos[imo] : -1;
301 Int_t trackIt = (ks == 1) && fTrackIt;
93831bcf 302 PushTrack(trackIt, iparent, kf,
303 p[0], p[1], p[2], p[3],
304 origin[0], origin[1], origin[2],
305 tof,
306 polar[0], polar[1], polar[2],
e2054d85 307 kPPrimary, nt, fHerwig->GetEVWGT(), ks);
618e1dc0 308 KeepTrack(nt);
309 newPos[i]=nt;
310 } // end of if: selection of particle
311 } // end of for: particle loop
312 if (newPos) delete[] newPos;
618e1dc0 313 // MakeHeader();
618e1dc0 314 if (nc > 0) {
315 jev+=nc;
316 if (jev >= fNpart || fNpart == -1) {
317 fKineBias=Float_t(fNpart)/Float_t(fTrials);
618e1dc0 318 break;
319 }
320 }
321 }
322 SetHighWaterMark(nt);
323// adjust weight due to kinematic selection
324 AdjustWeights();
325// get cross-section
326 fXsection=fHerwig->GetAVWGT();
327}
328
329void AliGenHerwig::AdjustWeights()
330{
331// Adjust the weights after generation of all events
332 TParticle *part;
5d12ce38 333 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
618e1dc0 334 for (Int_t i=0; i<ntrack; i++) {
5d12ce38 335 part= gAlice->GetMCApp()->Particle(i);
618e1dc0 336 part->SetWeight(part->GetWeight()*fKineBias);
337 }
338}
339
340
341void AliGenHerwig::KeepFullEvent()
342{
343 fKeep=1;
344}
345
346Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
347{
348//
349// Looks recursively if one of the daughters has been selected
350//
351// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
352 Int_t imin=-1;
353 Int_t imax=-1;
354 Int_t i;
355 Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
356 Bool_t selected=kFALSE;
357 if (hasDaughters) {
358 imin=iparticle->GetFirstDaughter();
359 imax=iparticle->GetLastDaughter();
360 for (i=imin; i<= imax; i++){
361 TParticle * jparticle = (TParticle *) particles->At(i);
362 Int_t ip=jparticle->GetPdgCode();
363 if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
364 selected=kTRUE; break;
365 }
366 if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
367 }
368 } else {
369 return kFALSE;
370 }
371
372 return selected;
373}
374
375
376Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
377{
378// Select flavor of particle
379// 0: all
380// 4: charm and beauty
381// 5: beauty
382 if (fFlavor == 0) return kTRUE;
383
384 Int_t ifl=TMath::Abs(pid/100);
385 if (ifl > 10) ifl/=10;
386 return (fFlavor == ifl);
387}
388
389Bool_t AliGenHerwig::Stable(TParticle* particle)
390{
391// Return true for a stable particle
392//
393 Int_t kf = TMath::Abs(particle->GetPdgCode());
394
395 if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
396
397 {
398 return kTRUE;
399 } else {
400 return kFALSE;
401 }
402}
403
404void AliGenHerwig::FinishRun()
405{
406 fHerwig->Hwefin();
407}
408
409
410AliGenHerwig& AliGenHerwig::operator=(const AliGenHerwig& rhs)
411{
412// Assignment operator
014a9521 413 rhs.Copy(*this);
414 return (*this);
618e1dc0 415}
416
417