ptharmax added.
[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.),
9305fc24 57 fPtHardMax(9999.),
2d654757 58 fPtRMS(0.),
a1d47d99 59 fMaxPr(10),
60 fMaxErrors(1000),
e2054d85 61 fEnSoft(1),
62 fEv1Pr(0),
2d654757 63 fEv2Pr(0),
9305fc24 64 fFileName(0),
65 fEtaMinParton(-20.),
66 fEtaMaxParton(20.),
67 fPhiMinParton(0.),
68 fPhiMaxParton(2.* TMath::Pi()),
69 fEtaMinGamma(-20.),
70 fEtaMaxGamma(20.),
71 fPhiMinGamma(0.),
72 fPhiMaxGamma(2. * TMath::Pi())
618e1dc0 73{
74// Constructor
fc7e1b1c 75 fEnergyCMS = 14000;
618e1dc0 76}
77
78AliGenHerwig::AliGenHerwig(Int_t npart)
2d654757 79 :AliGenMC(npart),
7677b281 80 fAutPDF("LHAPDF"),
81 fModPDF(19070),
2d654757 82 fStrucFunc(kCTEQ5L),
83 fKeep(0),
84 fDecaysOff(1),
85 fTrigger(0),
86 fSelectAll(0),
87 fFlavor(0),
2d654757 88 fMomentum1(7000),
89 fMomentum2(7000),
90 fKineBias(1),
91 fTrials(0),
92 fXsection(0),
93 fHerwig(0x0),
94 fProcess(0),
95 fPtHardMin(0.),
9305fc24 96 fPtHardMax(9999.),
2d654757 97 fPtRMS(0.),
98 fMaxPr(10),
99 fMaxErrors(1000),
100 fEnSoft(1),
101 fEv1Pr(0),
102 fEv2Pr(0),
9305fc24 103 fFileName(0),
104 fEtaMinParton(-20.),
105 fEtaMaxParton(20.),
106 fPhiMinParton(0.),
107 fPhiMaxParton(2.* TMath::Pi()),
108 fEtaMinGamma(-20.),
109 fEtaMaxGamma(20.),
110 fPhiMinGamma(0.),
111 fPhiMaxGamma(2. * TMath::Pi())
618e1dc0 112{
fc7e1b1c 113 fEnergyCMS = 14000;
618e1dc0 114 SetTarget();
115 SetProjectile();
7677b281 116 // Set random number generator
93831bcf 117 AliHerwigRndm::SetHerwigRandom(GetRandom());
618e1dc0 118}
119
618e1dc0 120AliGenHerwig::~AliGenHerwig()
121{
122// Destructor
123}
124
e2054d85 125void AliGenHerwig::SetEventListRange(Int_t eventFirst, Int_t eventLast)
126{
127 fEv1Pr = ++eventFirst;
128 fEv2Pr = ++eventLast;
129 if ( fEv2Pr == -1 ) fEv2Pr = fEv2Pr;
130}
131
618e1dc0 132void AliGenHerwig::Init()
133{
134// Initialisation
135 fTarget.Resize(8);
136 fProjectile.Resize(8);
137 SetMC(new THerwig6());
d274cd60 138 fHerwig=(THerwig6*) fMCEvGen;
618e1dc0 139 // initialize common blocks
140 fHerwig->Initialize(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
141 // reset parameters according to user needs
142 InitPDF();
143 fHerwig->SetPTMIN(fPtHardMin);
7a5495ad 144 fHerwig->SetPTMAX(fPtHardMax);
618e1dc0 145 fHerwig->SetPTRMS(fPtRMS);
146 fHerwig->SetMAXPR(fMaxPr);
147 fHerwig->SetMAXER(fMaxErrors);
148 fHerwig->SetENSOF(fEnSoft);
7677b281 149
e2054d85 150 fHerwig->SetEV1PR(fEv1Pr);
151 fHerwig->SetEV2PR(fEv2Pr);
152
153// C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
154// RMASS(1)=0.32
155// RMASS(2)=0.32
156// RMASS(3)=0.5
157// RMASS(4)=1.55
7677b281 158// RMASS(5)=4.75
e2054d85 159// RMASS(6)=174.3
160// RMASS(13)=0.75
161
162 fHerwig->SetRMASS(4,1.2);
163 fHerwig->SetRMASS(5,4.75);
7677b281 164
e2054d85 165 if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
166
9305fc24 167 //fHerwig->Hwusta("PI0 ");
e2054d85 168
618e1dc0 169 // compute parameter dependent constants
170 fHerwig->PrepareRun();
171}
172
7677b281 173void AliGenHerwig::InitJimmy()
174{
175// Initialisation
176 fTarget.Resize(8);
177 fProjectile.Resize(8);
178 SetMC(new THerwig6());
179 fHerwig=(THerwig6*) fMCEvGen;
180 // initialize common blocks
181 fHerwig->InitializeJimmy(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
182 // reset parameters according to user needs
183 InitPDF();
184 fHerwig->SetPTMIN(fPtHardMin);
185 fHerwig->SetPTRMS(fPtRMS);
186 fHerwig->SetMAXPR(fMaxPr);
187 fHerwig->SetMAXER(fMaxErrors);
188 fHerwig->SetENSOF(fEnSoft);
189
190 fHerwig->SetEV1PR(fEv1Pr);
191 fHerwig->SetEV2PR(fEv2Pr);
192
193// C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
194// RMASS(1)=0.32
195// RMASS(2)=0.32
196// RMASS(3)=0.5
197// RMASS(4)=1.55
198// RMASS(5)=4.75
199// RMASS(6)=174.3
200// RMASS(13)=0.75
201
202 fHerwig->SetRMASS(4,1.2);
203 fHerwig->SetRMASS(5,4.75);
204
205 if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
206
9305fc24 207 // fHerwig->Hwusta("PI0 ");
7677b281 208
209 // compute parameter dependent constants
210 fHerwig->PrepareRunJimmy();
211}
212
618e1dc0 213void AliGenHerwig::InitPDF()
214{
215 switch(fStrucFunc)
216 {
7677b281 217// ONLY USES LHAPDF STRUCTURE FUNCTIONS
618e1dc0 218 case kGRVLO98:
7677b281 219 fModPDF=80060;
220 fAutPDF="HWLHAPDF";
618e1dc0 221 break;
7677b281 222 case kCTEQ6:
223 fModPDF=10040;
224 fAutPDF="HWLHAPDF";
225 break;
226 case kCTEQ61:
227 fModPDF=10100;
228 fAutPDF="HWLHAPDF";
229 break;
230 case kCTEQ6m:
231 fModPDF=10050;
232 fAutPDF="HWLHAPDF";
233 break;
234 case kCTEQ6l:
235 fModPDF=10041;
236 fAutPDF="HWLHAPDF";
237 break;
238 case kCTEQ6ll:
239 fModPDF=10042;
240 fAutPDF="HWLHAPDF";
241 break;
242 case kCTEQ5M:
243 fModPDF=19050;
244 fAutPDF="HWLHAPDF";
618e1dc0 245 break;
246 case kCTEQ5L:
7677b281 247 fModPDF=19070;
248 fAutPDF="HWLHAPDF";
249 break;
250 case kCTEQ4M:
251 fModPDF=19150;
252 fAutPDF="HWLHAPDF";
253 break;
254 case kCTEQ4L:
255 fModPDF=19170;
256 fAutPDF="HWLHAPDF";
257 break;
56174e25 258// case kMRST2004nlo:
259// fModPDF=20400;
260// fAutPDF="HWLHAPDF";
261// break;
618e1dc0 262 default:
263 cerr << "This structure function is not inplemented " << fStrucFunc << endl;
264 break;
265 }
7677b281 266 fAutPDF.Resize(20);
618e1dc0 267 fHerwig->SetMODPDF(1,fModPDF);
268 fHerwig->SetMODPDF(2,fModPDF);
269 fHerwig->SetAUTPDF(1,fAutPDF);
270 fHerwig->SetAUTPDF(2,fAutPDF);
271}
272
273void AliGenHerwig::Generate()
274{
275 // Generate one event
276
277 Float_t polar[3] = {0,0,0};
278 Float_t origin[3]= {0,0,0};
279 Float_t origin0[3]= {0,0,0};
280 Float_t p[4], random[6];
7677b281 281
618e1dc0 282 static TClonesArray *particles;
283 // converts from mm/c to s
284 const Float_t kconv=0.001/2.999792458e8;
285 //
286 Int_t nt=0;
287 Int_t jev=0;
288 Int_t j, kf, ks, imo;
289 kf=0;
7677b281 290
618e1dc0 291 if(!particles) particles=new TClonesArray("TParticle",10000);
7677b281 292
618e1dc0 293 fTrials=0;
294 for (j=0;j<3;j++) origin0[j]=fOrigin[j];
295 if(fVertexSmear==kPerEvent) {
296 Rndm(random,6);
297 for (j=0;j<3;j++) {
298 origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
299 TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
300 }
301 }
7677b281 302
618e1dc0 303 while(1)
304 {
305 fHerwig->GenerateEvent();
306 fTrials++;
307 fHerwig->ImportParticles(particles,"All");
308 Int_t np = particles->GetEntriesFast()-1;
309 if (np == 0 ) continue;
7677b281 310
9305fc24 311 //Check hard partons or direct gamma in kine range
312
313 if (fProcess == kHeJets || fProcess == kHeDirectGamma) {
314 TParticle* parton1 = (TParticle *) particles->At(6);
315 TParticle* parton2 = (TParticle *) particles->At(7);
316 if (!CheckParton(parton1, parton2)) continue ;
317 }
318
618e1dc0 319 Int_t nc=0;
7677b281 320
618e1dc0 321 Int_t * newPos = new Int_t[np];
322 for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
7677b281 323
618e1dc0 324 for (Int_t i = 0; i<np; i++) {
325 TParticle * iparticle = (TParticle *) particles->At(i);
326 imo = iparticle->GetFirstMother();
327 kf = iparticle->GetPdgCode();
328 ks = iparticle->GetStatusCode();
7677b281 329 if (ks != 3 &&
330 KinematicSelection(iparticle,0))
618e1dc0 331 {
332 nc++;
333 p[0]=iparticle->Px();
334 p[1]=iparticle->Py();
335 p[2]=iparticle->Pz();
336 p[3]=iparticle->Energy();
337 origin[0]=origin0[0]+iparticle->Vx()/10;
338 origin[1]=origin0[1]+iparticle->Vy()/10;
339 origin[2]=origin0[2]+iparticle->Vz()/10;
340 Float_t tof = kconv*iparticle->T();
341 Int_t iparent = (imo > -1) ? newPos[imo] : -1;
342 Int_t trackIt = (ks == 1) && fTrackIt;
93831bcf 343 PushTrack(trackIt, iparent, kf,
344 p[0], p[1], p[2], p[3],
7677b281 345 origin[0], origin[1], origin[2],
93831bcf 346 tof,
347 polar[0], polar[1], polar[2],
e2054d85 348 kPPrimary, nt, fHerwig->GetEVWGT(), ks);
618e1dc0 349 KeepTrack(nt);
350 newPos[i]=nt;
351 } // end of if: selection of particle
352 } // end of for: particle loop
353 if (newPos) delete[] newPos;
618e1dc0 354 // MakeHeader();
618e1dc0 355 if (nc > 0) {
356 jev+=nc;
357 if (jev >= fNpart || fNpart == -1) {
358 fKineBias=Float_t(fNpart)/Float_t(fTrials);
618e1dc0 359 break;
360 }
361 }
362 }
363 SetHighWaterMark(nt);
364// adjust weight due to kinematic selection
365 AdjustWeights();
366// get cross-section
367 fXsection=fHerwig->GetAVWGT();
9305fc24 368 //printf(">> trials << %d\n",fTrials);
369}
370
371Bool_t AliGenHerwig::CheckParton(TParticle* parton1, TParticle* parton2)
372{
373// Check the kinematic trigger condition
374//
375//Select events with parton max energy
376 if(fPtHardMax < parton1->Pt()) return kFALSE;
377
378// Select events within angular window
379 Double_t eta[2];
380 eta[0] = parton1->Eta();
381 eta[1] = parton2->Eta();
382 Double_t phi[2];
383 phi[0] = parton1->Phi();
384 phi[1] = parton2->Phi();
385 Int_t pdg[2];
386 pdg[0] = parton1->GetPdgCode();
387 pdg[1] = parton2->GetPdgCode();
388 printf("min %f, max %f\n",fPtHardMin, fPtHardMax);
389 printf("Parton 1: %s, pT= %2.2f, eta = %1.2f, phi = %2.2f\n", parton1->GetName(),parton1->Pt(), eta[0],phi[0]*TMath::RadToDeg());
390 printf("Parton 2: %s, pT= %2.2f, eta = %1.2f, phi = %2.2f\n", parton2->GetName(),parton2->Pt(), eta[1],phi[1]*TMath::RadToDeg());
391
392 if (fProcess == kHeJets) {
393 //Check if one of the 2 outgoing partons are in the eta-phi window
394 for(Int_t i = 0; i < 2; i++)
395 if ((eta[i] < fEtaMaxParton && eta[i] > fEtaMinParton) &&
396 (phi[i] < fPhiMaxParton && phi[i] > fPhiMinParton)) return kTRUE ;
397 }
398
399 else {
400 //Check if the gamma and the jet are in the eta-phi window
401 Int_t igj = 0;
402 Int_t ijj = 0;
403 if(pdg[0] == 22) ijj=1;
404 else igj=1;
405 if ((eta[ijj] < fEtaMaxParton && eta[ijj] > fEtaMinParton) &&
406 (phi[ijj] < fPhiMaxParton && phi[ijj] > fPhiMinParton)) {
407
408 if ((eta[igj] < fEtaMaxGamma && eta[igj] > fEtaMinGamma) &&
409 (phi[igj] < fPhiMaxGamma && phi[igj] > fPhiMinGamma)) return kTRUE;
410
411 }
412 }
413
414 return kFALSE ;
618e1dc0 415}
416
417void AliGenHerwig::AdjustWeights()
418{
419// Adjust the weights after generation of all events
420 TParticle *part;
5d12ce38 421 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
618e1dc0 422 for (Int_t i=0; i<ntrack; i++) {
5d12ce38 423 part= gAlice->GetMCApp()->Particle(i);
618e1dc0 424 part->SetWeight(part->GetWeight()*fKineBias);
425 }
426}
7677b281 427
618e1dc0 428
429void AliGenHerwig::KeepFullEvent()
430{
431 fKeep=1;
432}
433
434Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
435{
436//
437// Looks recursively if one of the daughters has been selected
438//
439// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
440 Int_t imin=-1;
441 Int_t imax=-1;
442 Int_t i;
443 Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
444 Bool_t selected=kFALSE;
445 if (hasDaughters) {
446 imin=iparticle->GetFirstDaughter();
7677b281 447 imax=iparticle->GetLastDaughter();
618e1dc0 448 for (i=imin; i<= imax; i++){
7677b281 449 TParticle * jparticle = (TParticle *) particles->At(i);
618e1dc0 450 Int_t ip=jparticle->GetPdgCode();
451 if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
452 selected=kTRUE; break;
453 }
454 if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
455 }
456 } else {
457 return kFALSE;
458 }
459
460 return selected;
461}
462
463
464Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
465{
466// Select flavor of particle
467// 0: all
468// 4: charm and beauty
469// 5: beauty
470 if (fFlavor == 0) return kTRUE;
7677b281 471
618e1dc0 472 Int_t ifl=TMath::Abs(pid/100);
473 if (ifl > 10) ifl/=10;
474 return (fFlavor == ifl);
475}
476
477Bool_t AliGenHerwig::Stable(TParticle* particle)
478{
479// Return true for a stable particle
480//
481 Int_t kf = TMath::Abs(particle->GetPdgCode());
7677b281 482
618e1dc0 483 if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
7677b281 484
618e1dc0 485 {
486 return kTRUE;
487 } else {
488 return kFALSE;
489 }
490}
491
492void AliGenHerwig::FinishRun()
493{
494 fHerwig->Hwefin();
495}
496
7677b281 497void AliGenHerwig::FinishRunJimmy()
498{
499 fHerwig->Hwefin();
500 fHerwig->Jmefin();
501
502}
618e1dc0 503