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