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