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