Adding AliSelector, AliSelectorRL to ESD, STEER that can be used as base classes
[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
101AliGenHerwig::AliGenHerwig(const AliGenHerwig & Herwig)
2d654757 102 :AliGenMC(Herwig),
7677b281 103 fAutPDF("LHAPDF"),
104 fModPDF(19070),
2d654757 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);
7677b281 161
e2054d85 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
7677b281 170// RMASS(5)=4.75
e2054d85 171// RMASS(6)=174.3
172// RMASS(13)=0.75
173
174 fHerwig->SetRMASS(4,1.2);
175 fHerwig->SetRMASS(5,4.75);
7677b281 176
e2054d85 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
7677b281 185void AliGenHerwig::InitJimmy()
186{
187// Initialisation
188 fTarget.Resize(8);
189 fProjectile.Resize(8);
190 SetMC(new THerwig6());
191 fHerwig=(THerwig6*) fMCEvGen;
192 // initialize common blocks
193 fHerwig->InitializeJimmy(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
194 // reset parameters according to user needs
195 InitPDF();
196 fHerwig->SetPTMIN(fPtHardMin);
197 fHerwig->SetPTRMS(fPtRMS);
198 fHerwig->SetMAXPR(fMaxPr);
199 fHerwig->SetMAXER(fMaxErrors);
200 fHerwig->SetENSOF(fEnSoft);
201
202 fHerwig->SetEV1PR(fEv1Pr);
203 fHerwig->SetEV2PR(fEv2Pr);
204
205// C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
206// RMASS(1)=0.32
207// RMASS(2)=0.32
208// RMASS(3)=0.5
209// RMASS(4)=1.55
210// RMASS(5)=4.75
211// RMASS(6)=174.3
212// RMASS(13)=0.75
213
214 fHerwig->SetRMASS(4,1.2);
215 fHerwig->SetRMASS(5,4.75);
216
217 if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
218
219 fHerwig->Hwusta("PI0 ");
220
221 // compute parameter dependent constants
222 fHerwig->PrepareRunJimmy();
223}
224
618e1dc0 225void AliGenHerwig::InitPDF()
226{
227 switch(fStrucFunc)
228 {
7677b281 229// ONLY USES LHAPDF STRUCTURE FUNCTIONS
618e1dc0 230 case kGRVLO98:
7677b281 231 fModPDF=80060;
232 fAutPDF="HWLHAPDF";
618e1dc0 233 break;
7677b281 234 case kCTEQ6:
235 fModPDF=10040;
236 fAutPDF="HWLHAPDF";
237 break;
238 case kCTEQ61:
239 fModPDF=10100;
240 fAutPDF="HWLHAPDF";
241 break;
242 case kCTEQ6m:
243 fModPDF=10050;
244 fAutPDF="HWLHAPDF";
245 break;
246 case kCTEQ6l:
247 fModPDF=10041;
248 fAutPDF="HWLHAPDF";
249 break;
250 case kCTEQ6ll:
251 fModPDF=10042;
252 fAutPDF="HWLHAPDF";
253 break;
254 case kCTEQ5M:
255 fModPDF=19050;
256 fAutPDF="HWLHAPDF";
618e1dc0 257 break;
258 case kCTEQ5L:
7677b281 259 fModPDF=19070;
260 fAutPDF="HWLHAPDF";
261 break;
262 case kCTEQ4M:
263 fModPDF=19150;
264 fAutPDF="HWLHAPDF";
265 break;
266 case kCTEQ4L:
267 fModPDF=19170;
268 fAutPDF="HWLHAPDF";
269 break;
270 case kMRST2004nlo:
271 fModPDF=20400;
272 fAutPDF="HWLHAPDF";
618e1dc0 273 break;
274 default:
275 cerr << "This structure function is not inplemented " << fStrucFunc << endl;
276 break;
277 }
7677b281 278 fAutPDF.Resize(20);
618e1dc0 279 fHerwig->SetMODPDF(1,fModPDF);
280 fHerwig->SetMODPDF(2,fModPDF);
281 fHerwig->SetAUTPDF(1,fAutPDF);
282 fHerwig->SetAUTPDF(2,fAutPDF);
283}
284
285void AliGenHerwig::Generate()
286{
287 // Generate one event
288
289 Float_t polar[3] = {0,0,0};
290 Float_t origin[3]= {0,0,0};
291 Float_t origin0[3]= {0,0,0};
292 Float_t p[4], random[6];
7677b281 293
618e1dc0 294 static TClonesArray *particles;
295 // converts from mm/c to s
296 const Float_t kconv=0.001/2.999792458e8;
297 //
298 Int_t nt=0;
299 Int_t jev=0;
300 Int_t j, kf, ks, imo;
301 kf=0;
7677b281 302
618e1dc0 303 if(!particles) particles=new TClonesArray("TParticle",10000);
7677b281 304
618e1dc0 305 fTrials=0;
306 for (j=0;j<3;j++) origin0[j]=fOrigin[j];
307 if(fVertexSmear==kPerEvent) {
308 Rndm(random,6);
309 for (j=0;j<3;j++) {
310 origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
311 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
312 }
313 }
7677b281 314
618e1dc0 315 while(1)
316 {
317 fHerwig->GenerateEvent();
318 fTrials++;
319 fHerwig->ImportParticles(particles,"All");
320 Int_t np = particles->GetEntriesFast()-1;
321 if (np == 0 ) continue;
7677b281 322
618e1dc0 323 Int_t nc=0;
7677b281 324
618e1dc0 325 Int_t * newPos = new Int_t[np];
326 for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
7677b281 327
618e1dc0 328 for (Int_t i = 0; i<np; i++) {
329 TParticle * iparticle = (TParticle *) particles->At(i);
330 imo = iparticle->GetFirstMother();
331 kf = iparticle->GetPdgCode();
332 ks = iparticle->GetStatusCode();
7677b281 333 if (ks != 3 &&
334 KinematicSelection(iparticle,0))
618e1dc0 335 {
336 nc++;
337 p[0]=iparticle->Px();
338 p[1]=iparticle->Py();
339 p[2]=iparticle->Pz();
340 p[3]=iparticle->Energy();
341 origin[0]=origin0[0]+iparticle->Vx()/10;
342 origin[1]=origin0[1]+iparticle->Vy()/10;
343 origin[2]=origin0[2]+iparticle->Vz()/10;
344 Float_t tof = kconv*iparticle->T();
345 Int_t iparent = (imo > -1) ? newPos[imo] : -1;
346 Int_t trackIt = (ks == 1) && fTrackIt;
93831bcf 347 PushTrack(trackIt, iparent, kf,
348 p[0], p[1], p[2], p[3],
7677b281 349 origin[0], origin[1], origin[2],
93831bcf 350 tof,
351 polar[0], polar[1], polar[2],
e2054d85 352 kPPrimary, nt, fHerwig->GetEVWGT(), ks);
618e1dc0 353 KeepTrack(nt);
354 newPos[i]=nt;
355 } // end of if: selection of particle
356 } // end of for: particle loop
357 if (newPos) delete[] newPos;
618e1dc0 358 // MakeHeader();
618e1dc0 359 if (nc > 0) {
360 jev+=nc;
361 if (jev >= fNpart || fNpart == -1) {
362 fKineBias=Float_t(fNpart)/Float_t(fTrials);
618e1dc0 363 break;
364 }
365 }
366 }
367 SetHighWaterMark(nt);
368// adjust weight due to kinematic selection
369 AdjustWeights();
370// get cross-section
371 fXsection=fHerwig->GetAVWGT();
372}
373
374void AliGenHerwig::AdjustWeights()
375{
376// Adjust the weights after generation of all events
377 TParticle *part;
5d12ce38 378 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
618e1dc0 379 for (Int_t i=0; i<ntrack; i++) {
5d12ce38 380 part= gAlice->GetMCApp()->Particle(i);
618e1dc0 381 part->SetWeight(part->GetWeight()*fKineBias);
382 }
383}
7677b281 384
618e1dc0 385
386void AliGenHerwig::KeepFullEvent()
387{
388 fKeep=1;
389}
390
391Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
392{
393//
394// Looks recursively if one of the daughters has been selected
395//
396// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
397 Int_t imin=-1;
398 Int_t imax=-1;
399 Int_t i;
400 Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
401 Bool_t selected=kFALSE;
402 if (hasDaughters) {
403 imin=iparticle->GetFirstDaughter();
7677b281 404 imax=iparticle->GetLastDaughter();
618e1dc0 405 for (i=imin; i<= imax; i++){
7677b281 406 TParticle * jparticle = (TParticle *) particles->At(i);
618e1dc0 407 Int_t ip=jparticle->GetPdgCode();
408 if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
409 selected=kTRUE; break;
410 }
411 if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
412 }
413 } else {
414 return kFALSE;
415 }
416
417 return selected;
418}
419
420
421Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
422{
423// Select flavor of particle
424// 0: all
425// 4: charm and beauty
426// 5: beauty
427 if (fFlavor == 0) return kTRUE;
7677b281 428
618e1dc0 429 Int_t ifl=TMath::Abs(pid/100);
430 if (ifl > 10) ifl/=10;
431 return (fFlavor == ifl);
432}
433
434Bool_t AliGenHerwig::Stable(TParticle* particle)
435{
436// Return true for a stable particle
437//
438 Int_t kf = TMath::Abs(particle->GetPdgCode());
7677b281 439
618e1dc0 440 if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
7677b281 441
618e1dc0 442 {
443 return kTRUE;
444 } else {
445 return kFALSE;
446 }
447}
448
449void AliGenHerwig::FinishRun()
450{
451 fHerwig->Hwefin();
452}
453
454
455AliGenHerwig& AliGenHerwig::operator=(const AliGenHerwig& rhs)
456{
457// Assignment operator
014a9521 458 rhs.Copy(*this);
459 return (*this);
618e1dc0 460}
461
7677b281 462void AliGenHerwig::FinishRunJimmy()
463{
464 fHerwig->Hwefin();
465 fHerwig->Jmefin();
466
467}
618e1dc0 468