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