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