Moving to the new VMC naming convention
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx
CommitLineData
8d2cd130 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$ */
8d2cd130 17
18//
19// Generator using the TPythia interface (via AliPythia)
20// to generate pp collisions.
21// Using SetNuclei() also nuclear modifications to the structure functions
22// can be taken into account. This makes, of course, only sense for the
23// generation of the products of hard processes (heavy flavor, jets ...)
24//
25// andreas.morsch@cern.ch
26//
27
28#include <TDatabasePDG.h>
29#include <TParticle.h>
30#include <TPDGCode.h>
31#include <TSystem.h>
32#include <TTree.h>
33
34#include "AliConst.h"
35#include "AliDecayerPythia.h"
36#include "AliGenPythia.h"
37#include "AliGenPythiaEventHeader.h"
38#include "AliPythia.h"
7cdba479 39#include "AliPythiaRndm.h"
8d2cd130 40#include "AliRun.h"
41
42 ClassImp(AliGenPythia)
43
44AliGenPythia::AliGenPythia()
45 :AliGenMC()
46{
47// Default Constructor
48 fParticles = 0;
49 fPythia = 0;
e5c87a3d 50 fHeader = 0;
8d2cd130 51 fDecayer = new AliDecayerPythia();
52 SetEventListRange();
53 SetJetPhiRange();
54 SetJetEtaRange();
55 SetJetEtRange();
56 SetGammaPhiRange();
57 SetGammaEtaRange();
58 SetPtKick();
1d568bc2 59 fSetNuclei = kFALSE;
7cdba479 60 if (!AliPythiaRndm::GetPythiaRandom())
61 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 62}
63
64AliGenPythia::AliGenPythia(Int_t npart)
65 :AliGenMC(npart)
66{
67// default charm production at 5. 5 TeV
68// semimuonic decay
69// structure function GRVHO
70//
71 fName = "Pythia";
72 fTitle= "Particle Generator using PYTHIA";
73 fXsection = 0.;
8d2cd130 74 SetProcess();
75 SetStrucFunc();
76 SetForceDecay();
77 SetPtHard();
78 SetYHard();
79 SetEnergyCMS();
80 fDecayer = new AliDecayerPythia();
81 // Set random number generator
7cdba479 82 if (!AliPythiaRndm::GetPythiaRandom())
83 AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130 84 fFlavorSelect = 0;
85 // Produced particles
86 fParticles = new TClonesArray("TParticle",1000);
e5c87a3d 87 fHeader = 0;
8d2cd130 88 SetEventListRange();
89 SetJetPhiRange();
90 SetJetEtaRange();
91 SetJetEtRange();
92 SetGammaPhiRange();
93 SetGammaEtaRange();
94 SetJetReconstructionMode();
95 SetPtKick();
96 // Options determining what to keep in the stack (Heavy flavour generation)
97 fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
98 fFeedDownOpt = kTRUE; // allow feed down from higher family
99 // Fragmentation on/off
100 fFragmentation = kTRUE;
101 // Default counting mode
102 fCountMode = kCountAll;
592f8307 103 // Pycel
104 SetPycellParameters();
1d568bc2 105 fSetNuclei = kFALSE;
8d2cd130 106}
107
108AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
109{
110// copy constructor
111 Pythia.Copy(*this);
112}
113
114AliGenPythia::~AliGenPythia()
115{
116// Destructor
117}
118
592f8307 119void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
120 Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
121{
122// Set pycell parameters
123 fPycellEtaMax = etamax;
124 fPycellNEta = neta;
125 fPycellNPhi = nphi;
126 fPycellThreshold = thresh;
127 fPycellEtSeed = etseed;
128 fPycellMinEtJet = minet;
129 fPycellMaxRadius = r;
130}
131
132
133
8d2cd130 134void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
135{
136 // Set a range of event numbers, for which a table
137 // of generated particle will be printed
138 fDebugEventFirst = eventFirst;
139 fDebugEventLast = eventLast;
140 if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
141}
142
143void AliGenPythia::Init()
144{
145// Initialisation
146
147 SetMC(AliPythia::Instance());
148 fPythia=(AliPythia*) fgMCEvGen;
149//
150 fParentWeight=1./Float_t(fNpart);
151//
152// Forward Paramters to the AliPythia object
153 fDecayer->SetForceDecay(fForceDecay);
154 fDecayer->Init();
155
156
157 fPythia->SetCKIN(3,fPtHardMin);
158 fPythia->SetCKIN(4,fPtHardMax);
159 fPythia->SetCKIN(7,fYHardMin);
160 fPythia->SetCKIN(8,fYHardMax);
161
1a626d4e 162 if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);
8d2cd130 163 // Fragmentation?
164 if (fFragmentation) {
165 fPythia->SetMSTP(111,1);
166 } else {
167 fPythia->SetMSTP(111,0);
168 }
169
170
171// initial state radiation
172 fPythia->SetMSTP(61,fGinit);
173// final state radiation
174 fPythia->SetMSTP(71,fGfinal);
175// pt - kick
176 if (fPtKick > 0.) {
177 fPythia->SetMSTP(91,1);
178 fPythia->SetPARP(91,fPtKick);
179 } else {
180 fPythia->SetMSTP(91,0);
181 }
182
183 // fPythia->SetMSTJ(1,2);
184 //
185 fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
186
187// Parent and Children Selection
188 switch (fProcess)
189 {
190 case kPyCharm:
191 case kPyCharmUnforced:
adf4d898 192 case kPyCharmPbPbMNR:
aabc7187 193 case kPyCharmppMNR:
194 case kPyCharmpPbMNR:
8d2cd130 195 fParentSelect[0] = 411;
196 fParentSelect[1] = 421;
197 fParentSelect[2] = 431;
198 fParentSelect[3] = 4122;
199 fFlavorSelect = 4;
200 break;
adf4d898 201 case kPyD0PbPbMNR:
202 case kPyD0pPbMNR:
203 case kPyD0ppMNR:
8d2cd130 204 fParentSelect[0] = 421;
205 fFlavorSelect = 4;
206 break;
207 case kPyBeauty:
adf4d898 208 case kPyBeautyPbPbMNR:
209 case kPyBeautypPbMNR:
210 case kPyBeautyppMNR:
8d2cd130 211 fParentSelect[0]= 511;
212 fParentSelect[1]= 521;
213 fParentSelect[2]= 531;
214 fParentSelect[3]= 5122;
215 fParentSelect[4]= 5132;
216 fParentSelect[5]= 5232;
217 fParentSelect[6]= 5332;
218 fFlavorSelect = 5;
219 break;
220 case kPyBeautyUnforced:
221 fParentSelect[0] = 511;
222 fParentSelect[1] = 521;
223 fParentSelect[2] = 531;
224 fParentSelect[3] = 5122;
225 fParentSelect[4] = 5132;
226 fParentSelect[5] = 5232;
227 fParentSelect[6] = 5332;
228 fFlavorSelect = 5;
229 break;
230 case kPyJpsiChi:
231 case kPyJpsi:
232 fParentSelect[0] = 443;
233 break;
234 case kPyMb:
235 case kPyMbNonDiffr:
236 case kPyJets:
237 case kPyDirectGamma:
238 break;
239 }
240//
592f8307 241//
242// JetFinder for Trigger
243//
244// Configure detector (EMCAL like)
245//
246 fPythia->SetPARU(51, fPycellEtaMax);
247 fPythia->SetMSTU(51, fPycellNEta);
248 fPythia->SetMSTU(52, fPycellNPhi);
249//
250// Configure Jet Finder
251//
252 fPythia->SetPARU(58, fPycellThreshold);
253 fPythia->SetPARU(52, fPycellEtSeed);
254 fPythia->SetPARU(53, fPycellMinEtJet);
255 fPythia->SetPARU(54, fPycellMaxRadius);
256 fPythia->SetMSTU(54, 2);
257//
8d2cd130 258// This counts the total number of calls to Pyevnt() per run.
259 fTrialsRun = 0;
260 fQ = 0.;
261 fX1 = 0.;
262 fX2 = 0.;
263 fNev = 0 ;
264//
1d568bc2 265//
266//
8d2cd130 267 AliGenMC::Init();
1d568bc2 268//
269//
270//
271 if (fSetNuclei) {
272 fDyBoost = 0;
273 Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
274 }
8d2cd130 275}
276
277void AliGenPythia::Generate()
278{
279// Generate one event
e5c87a3d 280
8d2cd130 281 fDecayer->ForceDecay();
282
283 Float_t polar[3] = {0,0,0};
284 Float_t origin[3] = {0,0,0};
285 Float_t p[3];
286// converts from mm/c to s
287 const Float_t kconv=0.001/2.999792458e8;
288//
289 Int_t nt=0;
290 Int_t jev=0;
291 Int_t j, kf;
292 fTrials=0;
293
294 // Set collision vertex position
295 if(fVertexSmear==kPerEvent) {
296 fPythia->SetMSTP(151,1);
297 for (j=0;j<3;j++) {
298 fPythia->SetPARP(151+j, fOsigma[j]*10.);
299 }
300 } else if (fVertexSmear==kPerTrack) {
301 fPythia->SetMSTP(151,0);
302 }
303// event loop
304 while(1)
305 {
306 fPythia->Pyevnt();
307 if (gAlice->GetEvNumber()>=fDebugEventFirst &&
308 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
309 fTrials++;
310
311 fPythia->ImportParticles(fParticles,"All");
312
1d568bc2 313 Boost();
8d2cd130 314//
315//
316//
317 Int_t i;
318
319 Int_t np = fParticles->GetEntriesFast();
320 if (np == 0 ) continue;
321// Get event vertex and discard the event if the Z coord. is too big
322 TParticle *iparticle = (TParticle *) fParticles->At(0);
323 Float_t distz = iparticle->Vz()/10.;
324 if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue;
325//
d25cfd65 326 fVertex[0] = iparticle->Vx()/10.+fOrigin.At(0);
327 fVertex[1] = iparticle->Vy()/10.+fOrigin.At(1);
328 fVertex[2] = iparticle->Vz()/10.+fOrigin.At(2);
8d2cd130 329//
330 Int_t* pParent = new Int_t[np];
331 Int_t* pSelected = new Int_t[np];
332 Int_t* trackIt = new Int_t[np];
333 for (i=0; i< np; i++) {
334 pParent[i] = -1;
335 pSelected[i] = 0;
336 trackIt[i] = 0;
337 }
338
339 Int_t nc = 0; // Total n. of selected particles
340 Int_t nParents = 0; // Selected parents
341 Int_t nTkbles = 0; // Trackable particles
342 if (fProcess != kPyMb && fProcess != kPyJets &&
343 fProcess != kPyDirectGamma &&
344 fProcess != kPyMbNonDiffr) {
345
346 for (i = 0; i<np; i++) {
347 iparticle = (TParticle *) fParticles->At(i);
348 Int_t ks = iparticle->GetStatusCode();
349 kf = CheckPDGCode(iparticle->GetPdgCode());
350// No initial state partons
351 if (ks==21) continue;
352//
353// Heavy Flavor Selection
354//
355 // quark ?
356 kf = TMath::Abs(kf);
357 Int_t kfl = kf;
358 // meson ?
359 if (kfl > 10) kfl/=100;
360 // baryon
361 if (kfl > 10) kfl/=10;
362 if (kfl > 10) kfl/=10;
363
364 Int_t ipa = iparticle->GetFirstMother()-1;
365 Int_t kfMo = 0;
366
367 if (ipa > -1) {
368 TParticle * mother = (TParticle *) fParticles->At(ipa);
369 kfMo = TMath::Abs(mother->GetPdgCode());
370 }
371 // What to keep in Stack?
372 Bool_t flavorOK = kFALSE;
373 Bool_t selectOK = kFALSE;
374 if (fFeedDownOpt) {
375 if (kfl >= fFlavorSelect) flavorOK = kTRUE;
376 } else {
377 if (kfl > fFlavorSelect) {
378 nc = -1;
379 break;
380 }
381 if (kfl == fFlavorSelect) flavorOK = kTRUE;
382 }
383 switch (fStackFillOpt) {
384 case kFlavorSelection:
385 selectOK = kTRUE;
386 break;
387 case kParentSelection:
388 if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
389 break;
390 }
391 if (flavorOK && selectOK) {
392//
393// Heavy flavor hadron or quark
394//
395// Kinematic seletion on final state heavy flavor mesons
396 if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
397 {
398 continue;
399 }
400 pSelected[i] = 1;
401 if (ParentSelected(kf)) ++nParents; // Update parent count
402// printf("\n particle (HF) %d %d %d", i, pSelected[i], kf);
403 } else {
404// Kinematic seletion on decay products
405 if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
406 && !KinematicSelection(iparticle, 1))
407 {
408 continue;
409 }
410//
411// Decay products
412// Select if mother was selected and is not tracked
413
414 if (pSelected[ipa] &&
415 !trackIt[ipa] && // mother will be tracked ?
416 kfMo != 5 && // mother is b-quark, don't store fragments
417 kfMo != 4 && // mother is c-quark, don't store fragments
418 kf != 92) // don't store string
419 {
420//
421// Semi-stable or de-selected: diselect decay products:
422//
423//
424 if (pSelected[i] == -1 || fDecayer->GetLifetime(kf) > fMaxLifeTime)
425 {
426 Int_t ipF = iparticle->GetFirstDaughter();
427 Int_t ipL = iparticle->GetLastDaughter();
428 if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
429 }
430// printf("\n particle (decay) %d %d %d", i, pSelected[i], kf);
431 pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
432 }
433 }
434 if (pSelected[i] == -1) pSelected[i] = 0;
435 if (!pSelected[i]) continue;
436 // Count quarks only if you did not include fragmentation
437 if (fFragmentation && kf <= 10) continue;
438 nc++;
439// Decision on tracking
440 trackIt[i] = 0;
441//
442// Track final state particle
443 if (ks == 1) trackIt[i] = 1;
444// Track semi-stable particles
d25cfd65 445 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime)) trackIt[i] = 1;
8d2cd130 446// Track particles selected by process if undecayed.
447 if (fForceDecay == kNoDecay) {
448 if (ParentSelected(kf)) trackIt[i] = 1;
449 } else {
450 if (ParentSelected(kf)) trackIt[i] = 0;
451 }
452 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
453//
454//
455
456 } // particle selection loop
457 if (nc > 0) {
458 for (i = 0; i<np; i++) {
459 if (!pSelected[i]) continue;
460 TParticle * iparticle = (TParticle *) fParticles->At(i);
461 kf = CheckPDGCode(iparticle->GetPdgCode());
462 Int_t ks = iparticle->GetStatusCode();
463 p[0] = iparticle->Px();
464 p[1] = iparticle->Py();
465 p[2] = iparticle->Pz();
466 origin[0] = fOrigin[0]+iparticle->Vx()/10.;
467 origin[1] = fOrigin[1]+iparticle->Vy()/10.;
468 origin[2] = fOrigin[2]+iparticle->Vz()/10.;
469 Float_t tof = kconv*iparticle->T();
470 Int_t ipa = iparticle->GetFirstMother()-1;
471 Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
642f15cf 472 PushTrack(fTrackIt*trackIt[i] ,
8d2cd130 473 iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1., ks);
474 pParent[i] = nt;
475 KeepTrack(nt);
642f15cf 476 } // PushTrack loop
8d2cd130 477 }
478 } else {
479 nc = GenerateMB();
480 } // mb ?
481
482 if (pParent) delete[] pParent;
483 if (pSelected) delete[] pSelected;
484 if (trackIt) delete[] trackIt;
485
486 if (nc > 0) {
487 switch (fCountMode) {
488 case kCountAll:
489 // printf(" Count all \n");
490 jev += nc;
491 break;
492 case kCountParents:
493 // printf(" Count parents \n");
494 jev += nParents;
495 break;
496 case kCountTrackables:
497 // printf(" Count trackable \n");
498 jev += nTkbles;
499 break;
500 }
501 if (jev >= fNpart || fNpart == -1) {
502 fKineBias=Float_t(fNpart)/Float_t(fTrials);
503 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
504
505 fQ += fPythia->GetVINT(51);
506 fX1 += fPythia->GetVINT(41);
507 fX2 += fPythia->GetVINT(42);
508 fTrialsRun += fTrials;
509 fNev++;
510 MakeHeader();
511 break;
512 }
513 }
514 } // event loop
515 SetHighWaterMark(nt);
516// adjust weight due to kinematic selection
517 AdjustWeights();
518// get cross-section
519 fXsection=fPythia->GetPARI(1);
520}
521
522Int_t AliGenPythia::GenerateMB()
523{
524//
525// Min Bias selection and other global selections
526//
527 Int_t i, kf, nt, iparent;
528 Int_t nc = 0;
529 Float_t p[3];
530 Float_t polar[3] = {0,0,0};
531 Float_t origin[3] = {0,0,0};
532// converts from mm/c to s
533 const Float_t kconv=0.001/2.999792458e8;
534
535 Int_t np = fParticles->GetEntriesFast();
536 Int_t* pParent = new Int_t[np];
537 for (i=0; i< np; i++) pParent[i] = -1;
538 if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
539 TParticle* jet1 = (TParticle *) fParticles->At(6);
540 TParticle* jet2 = (TParticle *) fParticles->At(7);
541 if (!CheckTrigger(jet1, jet2)) return 0;
542 }
543
544 for (i = 0; i<np; i++) {
545 Int_t trackIt = 0;
546 TParticle * iparticle = (TParticle *) fParticles->At(i);
547 kf = CheckPDGCode(iparticle->GetPdgCode());
548 Int_t ks = iparticle->GetStatusCode();
549 Int_t km = iparticle->GetFirstMother();
550 if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
551 (ks != 1) ||
552 (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
553 nc++;
554 if (ks == 1) trackIt = 1;
555 Int_t ipa = iparticle->GetFirstMother()-1;
556
557 iparent = (ipa > -1) ? pParent[ipa] : -1;
558
559//
560// store track information
561 p[0] = iparticle->Px();
562 p[1] = iparticle->Py();
563 p[2] = iparticle->Pz();
564 origin[0] = fOrigin[0]+iparticle->Vx()/10.;
565 origin[1] = fOrigin[1]+iparticle->Vy()/10.;
566 origin[2] = fOrigin[2]+iparticle->Vz()/10.;
567 Float_t tof=kconv*iparticle->T();
642f15cf 568 PushTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
8d2cd130 569 tof, kPPrimary, nt, 1., ks);
570 KeepTrack(nt);
571 pParent[i] = nt;
572 } // select particle
573 } // particle loop
574
575 if (pParent) delete[] pParent;
576
577 printf("\n I've put %i particles on the stack \n",nc);
578 return nc;
579}
580
581
582void AliGenPythia::FinishRun()
583{
584// Print x-section summary
585 fPythia->Pystat(1);
586 fQ /= fNev;
587 fX1 /= fNev;
588 fX2 /= fNev;
589 printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
590 printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
591
592
593}
594
595void AliGenPythia::AdjustWeights()
596{
597// Adjust the weights after generation of all events
598//
599 TParticle *part;
600 Int_t ntrack=gAlice->GetNtrack();
601 for (Int_t i=0; i<ntrack; i++) {
602 part= gAlice->Particle(i);
603 part->SetWeight(part->GetWeight()*fKineBias);
604 }
605}
606
607void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
608{
609// Treat protons as inside nuclei with mass numbers a1 and a2
1d568bc2 610
1a626d4e 611 fAProjectile = a1;
612 fATarget = a2;
1d568bc2 613 fSetNuclei = kTRUE;
8d2cd130 614}
615
616
617void AliGenPythia::MakeHeader()
618{
619// Builds the event header, to be called after each event
e5c87a3d 620 if (fHeader) delete fHeader;
621 fHeader = new AliGenPythiaEventHeader("Pythia");
8d2cd130 622//
623// Event type
e5c87a3d 624 ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
8d2cd130 625//
626// Number of trials
e5c87a3d 627 ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
8d2cd130 628//
629// Event Vertex
d25cfd65 630 fHeader->SetPrimaryVertex(fVertex);
8d2cd130 631//
632// Jets that have triggered
633 if (fProcess == kPyJets)
634 {
635 Int_t ntrig, njet;
636 Float_t jets[4][10];
637 GetJets(njet, ntrig, jets);
638
639 for (Int_t i = 0; i < ntrig; i++) {
e5c87a3d 640 ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
8d2cd130 641 jets[3][i]);
642 }
643 }
e5c87a3d 644 gAlice->SetGenEventHeader(fHeader);
8d2cd130 645}
646
647
648Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
649{
650// Check the kinematic trigger condition
651//
652 Double_t eta[2];
653 eta[0] = jet1->Eta();
654 eta[1] = jet2->Eta();
655 Double_t phi[2];
656 phi[0] = jet1->Phi();
657 phi[1] = jet2->Phi();
658 Int_t pdg[2];
659 pdg[0] = jet1->GetPdgCode();
660 pdg[1] = jet2->GetPdgCode();
661 Bool_t triggered = kFALSE;
662
663 if (fProcess == kPyJets) {
664 Int_t njets = 0;
665 Int_t ntrig = 0;
666 Float_t jets[4][10];
667//
668// Use Pythia clustering on parton level to determine jet axis
669//
670 GetJets(njets, ntrig, jets);
671
672 if (ntrig) triggered = kTRUE;
673//
674 } else {
675 Int_t ij = 0;
676 Int_t ig = 1;
677 if (pdg[0] == kGamma) {
678 ij = 1;
679 ig = 0;
680 }
681 //Check eta range first...
682 if ((eta[ij] < fEtaMaxJet && eta[ij] > fEtaMinJet) &&
683 (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
684 {
685 //Eta is okay, now check phi range
686 if ((phi[ij] < fPhiMaxJet && phi[ij] > fPhiMinJet) &&
687 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
688 {
689 triggered = kTRUE;
690 }
691 }
692 }
693 return triggered;
694}
695
696AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs)
697{
698// Assignment operator
699 return *this;
700}
701
702void AliGenPythia::LoadEvent()
703{
704//
705// Load event into Pythia Common Block
706//
707
708
709 Int_t npart = (Int_t) (gAlice->TreeK())->GetEntries();
710 (fPythia->GetPyjets())->N = npart;
711
712 for (Int_t part = 0; part < npart; part++) {
713 TParticle *MPart = gAlice->Particle(part);
714 Int_t kf = MPart->GetPdgCode();
715 Int_t ks = MPart->GetStatusCode();
716 Float_t px = MPart->Px();
717 Float_t py = MPart->Py();
718 Float_t pz = MPart->Pz();
719 Float_t e = MPart->Energy();
720 Float_t p = TMath::Sqrt(px * px + py * py + pz * pz);
721 Float_t m = TMath::Sqrt(e * e - p * p);
722
723
724 (fPythia->GetPyjets())->P[0][part] = px;
725 (fPythia->GetPyjets())->P[1][part] = py;
726 (fPythia->GetPyjets())->P[2][part] = pz;
727 (fPythia->GetPyjets())->P[3][part] = e;
728 (fPythia->GetPyjets())->P[4][part] = m;
729
730 (fPythia->GetPyjets())->K[1][part] = kf;
731 (fPythia->GetPyjets())->K[0][part] = ks;
732 }
733}
734
735void AliGenPythia::RecJetsUA1(Float_t eCellMin, Float_t eCellSeed, Float_t eMin, Float_t rMax,
736 Int_t& njets, Float_t jets [4][50])
737{
738//
739// Calls the Pythia jet finding algorithm to find jets in the current event
740//
741//
8d2cd130 742//
743// Save jets
744 Int_t n = fPythia->GetN();
745
746//
747// Run Jet Finder
748 fPythia->Pycell(njets);
749 Int_t i;
750 for (i = 0; i < njets; i++) {
751 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
752 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
753 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
754 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
755
756 jets[0][i] = px;
757 jets[1][i] = py;
758 jets[2][i] = pz;
759 jets[3][i] = e;
760 }
761}
762
763
764
765void AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
766{
767//
768// Calls the Pythia clustering algorithm to find jets in the current event
769//
770 Int_t n = fPythia->GetN();
771 nJets = 0;
772 nJetsTrig = 0;
773 if (fJetReconstruction == kCluster) {
774//
775// Configure cluster algorithm
776//
777 fPythia->SetPARU(43, 2.);
778 fPythia->SetMSTU(41, 1);
779//
780// Call cluster algorithm
781//
782 fPythia->Pyclus(nJets);
783//
784// Loading jets from common block
785//
786 } else {
592f8307 787
8d2cd130 788//
789// Run Jet Finder
790 fPythia->Pycell(nJets);
791 }
792
793 Int_t i;
794 for (i = 0; i < nJets; i++) {
795 Float_t px = (fPythia->GetPyjets())->P[0][n+i];
796 Float_t py = (fPythia->GetPyjets())->P[1][n+i];
797 Float_t pz = (fPythia->GetPyjets())->P[2][n+i];
798 Float_t e = (fPythia->GetPyjets())->P[3][n+i];
799 Float_t pt = TMath::Sqrt(px * px + py * py);
800 Float_t phi = TMath::ATan2(py,px);
801 Float_t theta = TMath::ATan2(pt,pz);
802 Float_t et = e * TMath::Sin(theta);
803 Float_t eta = -TMath::Log(TMath::Tan(theta / 2.));
804
805 if (
806 eta > fEtaMinJet && eta < fEtaMaxJet &&
675eb105 807 phi > fPhiMinJet && phi < fPhiMaxJet &&
8d2cd130 808 et > fEtMinJet && et < fEtMaxJet
809 )
810 {
811 jets[0][nJetsTrig] = px;
812 jets[1][nJetsTrig] = py;
813 jets[2][nJetsTrig] = pz;
814 jets[3][nJetsTrig] = e;
815 nJetsTrig++;
816
817 } else {
818// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
819 }
820 }
821}
822
823
824#ifdef never
825void AliGenPythia::Streamer(TBuffer &R__b)
826{
827 // Stream an object of class AliGenPythia.
828
829 if (R__b.IsReading()) {
830 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
831 AliGenerator::Streamer(R__b);
832 R__b >> (Int_t&)fProcess;
833 R__b >> (Int_t&)fStrucFunc;
834 R__b >> (Int_t&)fForceDecay;
835 R__b >> fEnergyCMS;
836 R__b >> fKineBias;
837 R__b >> fTrials;
838 fParentSelect.Streamer(R__b);
839 fChildSelect.Streamer(R__b);
840 R__b >> fXsection;
841// (AliPythia::Instance())->Streamer(R__b);
842 R__b >> fPtHardMin;
843 R__b >> fPtHardMax;
844// if (fDecayer) fDecayer->Streamer(R__b);
845 } else {
846 R__b.WriteVersion(AliGenPythia::IsA());
847 AliGenerator::Streamer(R__b);
848 R__b << (Int_t)fProcess;
849 R__b << (Int_t)fStrucFunc;
850 R__b << (Int_t)fForceDecay;
851 R__b << fEnergyCMS;
852 R__b << fKineBias;
853 R__b << fTrials;
854 fParentSelect.Streamer(R__b);
855 fChildSelect.Streamer(R__b);
856 R__b << fXsection;
857// R__b << fPythia;
858 R__b << fPtHardMin;
859 R__b << fPtHardMax;
860 // fDecayer->Streamer(R__b);
861 }
862}
863#endif
864