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