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