]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - THerwig/AliGenHerwig.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
updated
[u/mrichter/AliRoot.git] / THerwig / AliGenHerwig.cxx
... / ...
CommitLineData
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
16/* $Id$ */
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
24#include <Riostream.h>
25#include <TClonesArray.h>
26#include <TParticle.h>
27
28#include <THerwig6.h>
29
30#include "AliGenHerwig.h"
31#include "AliGenHerwigEventHeader.h"
32#include "AliHerwigRndm.h"
33#include "AliMC.h"
34#include "AliRun.h"
35#include "driver.h"
36
37using std::cerr;
38using std::endl;
39
40ClassImp(AliGenHerwig)
41
42
43 AliGenHerwig::AliGenHerwig() :
44 AliGenMC(),
45 fAutPDF("LHAPDF"),
46 fModPDF(19070),
47 fStrucFunc(kCTEQ5L),
48 fKeep(0),
49 fDecaysOff(1),
50 fTrigger(0),
51 fSelectAll(0),
52 fFlavor(0),
53 fMomentum1(7000),
54 fMomentum2(7000),
55 fKineBias(1),
56 fTrials(0),
57 fXsection(0),
58 fHerwig(0x0),
59 fProcess(0),
60 fPtHardMin(0.),
61 fPtHardMax(9999.),
62 fPtRMS(0.),
63 fMaxPr(10),
64 fMaxErrors(1000),
65 fEnSoft(1),
66 fEv1Pr(0),
67 fEv2Pr(0),
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.),
76 fPhiMaxGamma(2. * TMath::Pi()),
77 fHeader(0)
78{
79// Constructor
80 fEnergyCMS = 14000;
81}
82
83AliGenHerwig::AliGenHerwig(Int_t npart)
84 :AliGenMC(npart),
85 fAutPDF("LHAPDF"),
86 fModPDF(19070),
87 fStrucFunc(kCTEQ5L),
88 fKeep(0),
89 fDecaysOff(1),
90 fTrigger(0),
91 fSelectAll(0),
92 fFlavor(0),
93 fMomentum1(7000),
94 fMomentum2(7000),
95 fKineBias(1),
96 fTrials(0),
97 fXsection(0),
98 fHerwig(0x0),
99 fProcess(0),
100 fPtHardMin(0.),
101 fPtHardMax(9999.),
102 fPtRMS(0.),
103 fMaxPr(10),
104 fMaxErrors(1000),
105 fEnSoft(1),
106 fEv1Pr(0),
107 fEv2Pr(0),
108 fFileName(0),
109 fEtaMinParton(-20.),
110 fEtaMaxParton(20.),
111 fPhiMinParton(0.),
112 fPhiMaxParton(2.* TMath::Pi()),
113 fEtaMinGamma(-20.),
114 fEtaMaxGamma(20.),
115 fPhiMinGamma(0.),
116 fPhiMaxGamma(2. * TMath::Pi()),
117 fHeader(0)
118{
119// Constructor
120 fEnergyCMS = 14000;
121 SetTarget();
122 SetProjectile();
123 // Set random number generator
124 AliHerwigRndm::SetHerwigRandom(GetRandom());
125}
126
127AliGenHerwig::~AliGenHerwig()
128{
129// Destructor
130}
131
132void AliGenHerwig::SetEventListRange(Int_t eventFirst, Int_t eventLast)
133{
134 fEv1Pr = eventFirst;
135 fEv2Pr = eventLast;
136 if ( fEv2Pr == -1 ) fEv2Pr = fEv1Pr;
137}
138
139void AliGenHerwig::Init()
140{
141// Initialisation
142 fTarget.Resize(8);
143 fProjectile.Resize(8);
144 SetMC(new THerwig6());
145 fHerwig=(THerwig6*) fMCEvGen;
146 // initialize common blocks
147 fHerwig->Initialize(fProjectile.Data(), fTarget.Data(), fMomentum1, fMomentum2, fProcess); //
148 // reset parameters according to user needs
149 InitPDF();
150 fHerwig->SetPTMIN(fPtHardMin);
151 fHerwig->SetPTMAX(fPtHardMax);
152 fHerwig->SetPTRMS(fPtRMS);
153 fHerwig->SetMAXPR(fMaxPr);
154 fHerwig->SetMAXER(fMaxErrors);
155 fHerwig->SetENSOF(fEnSoft);
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
161// RMASS(5)=4.75
162// RMASS(6)=174.3
163// RMASS(13)=0.75
164
165 fHerwig->SetRMASS(4,1.2);
166 fHerwig->SetRMASS(5,4.75);
167
168 if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(), 49);
169
170 //fHerwig->Hwusta("PI0 ");
171
172 // compute parameter dependent constants
173 fHerwig->PrepareRun();
174}
175
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
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(), 49);
206
207 // fHerwig->Hwusta("PI0 ");
208
209 // compute parameter dependent constants
210 fHerwig->PrepareRunJimmy();
211}
212
213void AliGenHerwig::InitPDF()
214{
215// Initialize PDF
216 switch(fStrucFunc)
217 {
218// ONLY USES LHAPDF STRUCTURE FUNCTIONS
219 case kGRVLO98:
220 fModPDF=80060;
221 fAutPDF="HWLHAPDF";
222 break;
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";
246 break;
247 case kCTEQ5L:
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;
259// case kMRST2004nlo:
260// fModPDF=20400;
261// fAutPDF="HWLHAPDF";
262// break;
263 default:
264 cerr << "This structure function is not inplemented " << fStrucFunc << endl;
265 break;
266 }
267 fAutPDF.Resize(20);
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
278 Float_t polar[3] = {0,0,0};
279 Float_t origin[3] = {0,0,0};
280 Float_t p[4];
281
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 kf, ks, imo;
289 kf=0;
290
291 if(!particles) particles=new TClonesArray("TParticle",10000);
292
293 fTrials=0;
294
295 // Set collision vertex position
296 if (fVertexSmear == kPerEvent) Vertex();
297
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;
305
306 //Check hard partons or direct gamma in kine range
307
308 if (fProcess == kHeJets || fProcess == kHeDirectGamma) {
309 TParticle* parton1 = (TParticle *) particles->At(6);
310 TParticle* parton2 = (TParticle *) particles->At(7);
311 if (!CheckParton(parton1, parton2)) continue ;
312 }
313
314 //
315 if (gAlice) {
316 if (gAlice->GetEvNumber()>=fEv1Pr &&
317 gAlice->GetEvNumber()<=fEv2Pr) fHerwig->PrintEvt();
318
319 }
320
321
322 Int_t nc = 0;
323 fNprimaries = 0;
324
325 Int_t * newPos = new Int_t[np];
326 for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
327
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();
333 if (ks != 3 &&
334 KinematicSelection(iparticle,0))
335 {
336 nc++;
337 p[0]=iparticle->Px();
338 p[1]=iparticle->Py();
339 p[2]=iparticle->Pz();
340 p[3]=iparticle->Energy();
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
346 Float_t tof = fTime + kconv*iparticle->T();
347 Int_t iparent = (imo > -1) ? newPos[imo] : -1;
348 Int_t trackIt = (ks == 1) && fTrackIt;
349 PushTrack(trackIt, iparent, kf,
350 p[0], p[1], p[2], p[3],
351 origin[0], origin[1], origin[2],
352 tof,
353 polar[0], polar[1], polar[2],
354 kPPrimary, nt, fHerwig->GetEVWGT(), ks);
355 KeepTrack(nt);
356 newPos[i]=nt;
357 fNprimaries++;
358 } // end of if: selection of particle
359 } // end of for: particle loop
360 if (newPos) delete[] newPos;
361 if (nc > 0) {
362 jev+=nc;
363 if (jev >= fNpart || fNpart == -1) {
364 fKineBias=Float_t(fNpart)/Float_t(fTrials);
365 break;
366 }
367 }
368 }
369//
370 MakeHeader();
371//
372 SetHighWaterMark(nt);
373// adjust weight due to kinematic selection
374 AdjustWeights();
375// get cross-section
376 fXsection=fHerwig->GetAVWGT();
377 //printf(">> trials << %d\n",fTrials);
378}
379
380Bool_t AliGenHerwig::CheckParton(const TParticle* parton1, const TParticle* parton2)
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 ;
424}
425
426void AliGenHerwig::AdjustWeights()
427{
428// Adjust the weights after generation of all events
429 TParticle *part;
430 Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
431 for (Int_t i=0; i<ntrack; i++) {
432 part= gAlice->GetMCApp()->Particle(i);
433 part->SetWeight(part->GetWeight()*fKineBias);
434 }
435}
436
437
438void AliGenHerwig::KeepFullEvent()
439{
440 fKeep=1;
441}
442
443Bool_t AliGenHerwig::DaughtersSelection(const TParticle* iparticle, const TClonesArray* particles)
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();
456 imax=iparticle->GetLastDaughter();
457 for (i=imin; i<= imax; i++){
458 TParticle * jparticle = (TParticle *) particles->At(i);
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
473Bool_t AliGenHerwig::SelectFlavor(Int_t pid) const
474{
475// Select flavor of particle
476// 0: all
477// 4: charm and beauty
478// 5: beauty
479 if (fFlavor == 0) return kTRUE;
480
481 Int_t ifl=TMath::Abs(pid/100);
482 if (ifl > 10) ifl/=10;
483 return (fFlavor == ifl);
484}
485
486Bool_t AliGenHerwig::Stable(const TParticle* particle) const
487{
488// Return true for a stable particle
489//
490 Int_t kf = TMath::Abs(particle->GetPdgCode());
491
492 if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
493
494 {
495 return kTRUE;
496 } else {
497 return kFALSE;
498 }
499}
500
501void AliGenHerwig::FinishRun()
502{
503 fHerwig->Hwefin();
504}
505
506void AliGenHerwig::FinishRunJimmy()
507{
508 fHerwig->Hwefin();
509 fHerwig->Jmefin();
510
511}
512
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//
528// Event weight (cross section)
529 ((AliGenHerwigEventHeader*) fHeader)->SetWeight(fHerwig->GetEVWGT());
530//
531// Event Vertex
532 fHeader->SetPrimaryVertex(fVertex);
533 fHeader->SetInteractionTime(fTime);
534//
535// Number of primaries
536 fHeader->SetNProduced(fNprimaries);
537// Pass header
538//
539 AddHeader(fHeader);
540 fHeader = 0x0;
541}