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