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