PYTHIA6/AliGenPythia.cxx

   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
  16 /* $Id$ */
  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 #include "AliConst.h"
  34 #include "AliDecayerPythia.h"
  35 #include "AliGenPythia.h"
  36 #include "AliHeader.h"
  37 #include "AliGenPythiaEventHeader.h"
  38 #include "AliPythia.h"
  39 #include "AliPythiaRndm.h"
  40 #include "AliRun.h"
  41 #include "AliStack.h"
  42 #include "AliRunLoader.h"
  43 #include "AliMC.h"
  44 #include "pyquenCommon.h"
  45
  46 ClassImp(AliGenPythia)
  47
  48
  49 AliGenPythia::AliGenPythia():
  50     AliGenMC(),
  51     fProcess(kPyCharm),
  52     fStrucFunc(kCTEQ5L),
  53     fEnergyCMS(5500.),
  54     fKineBias(0.),
  55     fTrials(0),
  56     fTrialsRun(0),
  57     fQ(0.),
  58     fX1(0.),
  59     fX2(0.),
  60     fEventTime(0.),
  61     fInteractionRate(0.),
  62     fTimeWindow(0.),
  63     fCurSubEvent(0),
  64     fEventsTime(0),
  65     fNev(0),
  66     fFlavorSelect(0),
  67     fXsection(0.),
  68     fPythia(0),
  69     fPtHardMin(0.),
  70     fPtHardMax(1.e4),
  71     fYHardMin(-1.e10),
  72     fYHardMax(1.e10),
  73     fGinit(1),
  74     fGfinal(1),
  75     fHadronisation(1),
  76     fNpartons(0),
  77     fReadFromFile(0),
  78     fQuench(0),
  79     fPtKick(1.),
  80     fFullEvent(kTRUE),
  81     fDecayer(new AliDecayerPythia()),
  82     fDebugEventFirst(-1),
  83     fDebugEventLast(-1),
  84     fEtMinJet(0.),
  85     fEtMaxJet(1.e4),
  86     fEtaMinJet(-20.),
  87     fEtaMaxJet(20.),
  88     fPhiMinJet(0.),
  89     fPhiMaxJet(2.* TMath::Pi()),
  90     fJetReconstruction(kCell),
  91     fEtaMinGamma(-20.),
  92     fEtaMaxGamma(20.),
  93     fPhiMinGamma(0.),
  94     fPhiMaxGamma(2. * TMath::Pi()),
  95     fPycellEtaMax(2.),
  96     fPycellNEta(274),
  97     fPycellNPhi(432),
  98     fPycellThreshold(0.),
  99     fPycellEtSeed(4.),
 100     fPycellMinEtJet(10.),
 101     fPycellMaxRadius(1.),
 102     fStackFillOpt(kFlavorSelection),
 103     fFeedDownOpt(kTRUE),
 104     fFragmentation(kTRUE),
 105     fSetNuclei(kFALSE),
 106     fNewMIS(kFALSE),
 107     fHFoff(kFALSE),
 108     fTriggerParticle(0),
 109     fTriggerEta(0.9),
 110     fCountMode(kCountAll),
 111     fHeader(0),
 112     fRL(0),
 113     fFileName(0)
 114 {
 115 // Default Constructor
 116   SetNuclei(0,0);
 117   if (!AliPythiaRndm::GetPythiaRandom())
 118       AliPythiaRndm::SetPythiaRandom(GetRandom());
 119 }
 120
 121 AliGenPythia::AliGenPythia(Int_t npart)
 122     :AliGenMC(npart),
 123      fProcess(kPyCharm),
 124      fStrucFunc(kCTEQ5L),
 125      fEnergyCMS(5500.),
 126      fKineBias(0.),
 127      fTrials(0),
 128      fTrialsRun(0),
 129      fQ(0.),
 130      fX1(0.),
 131      fX2(0.),
 132      fEventTime(0.),
 133      fInteractionRate(0.),
 134      fTimeWindow(0.),
 135      fCurSubEvent(0),
 136      fEventsTime(0),
 137      fNev(0),
 138      fFlavorSelect(0),
 139      fXsection(0.),
 140      fPythia(0),
 141      fPtHardMin(0.),
 142      fPtHardMax(1.e4),
 143      fYHardMin(-1.e10),
 144      fYHardMax(1.e10),
 145      fGinit(kTRUE),
 146      fGfinal(kTRUE),
 147      fHadronisation(kTRUE),
 148      fNpartons(0),
 149      fReadFromFile(kFALSE),
 150      fQuench(kFALSE),
 151      fPtKick(1.),
 152      fFullEvent(kTRUE),
 153      fDecayer(new AliDecayerPythia()),
 154      fDebugEventFirst(-1),
 155      fDebugEventLast(-1),
 156      fEtMinJet(0.),
 157      fEtMaxJet(1.e4),
 158      fEtaMinJet(-20.),
 159      fEtaMaxJet(20.),
 160      fPhiMinJet(0.),
 161      fPhiMaxJet(2.* TMath::Pi()),
 162      fJetReconstruction(kCell),
 163      fEtaMinGamma(-20.),
 164      fEtaMaxGamma(20.),
 165      fPhiMinGamma(0.),
 166      fPhiMaxGamma(2. * TMath::Pi()),
 167      fPycellEtaMax(2.),
 168      fPycellNEta(274),
 169      fPycellNPhi(432),
 170      fPycellThreshold(0.),
 171      fPycellEtSeed(4.),
 172      fPycellMinEtJet(10.),
 173      fPycellMaxRadius(1.),
 174      fStackFillOpt(kFlavorSelection),
 175      fFeedDownOpt(kTRUE),
 176      fFragmentation(kTRUE),
 177      fSetNuclei(kFALSE),
 178      fNewMIS(kFALSE),
 179      fHFoff(kFALSE),
 180      fTriggerParticle(0),
 181      fTriggerEta(0.9),
 182      fCountMode(kCountAll),
 183      fHeader(0),
 184      fRL(0),
 185      fFileName(0)
 186 {
 187 // default charm production at 5. 5 TeV
 188 // semimuonic decay
 189 // structure function GRVHO
 190 //
 191     fName = "Pythia";
 192     fTitle= "Particle Generator using PYTHIA";
 193     SetForceDecay();
 194     // Set random number generator
 195     if (!AliPythiaRndm::GetPythiaRandom())
 196       AliPythiaRndm::SetPythiaRandom(GetRandom());
 197     fParticles = new TClonesArray("TParticle",1000);
 198     SetNuclei(0,0);
 199  }
 200
 201 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
 202     :AliGenMC(Pythia),
 203      fProcess(kPyCharm),
 204      fStrucFunc(kCTEQ5L),
 205      fEnergyCMS(5500.),
 206      fKineBias(0.),
 207      fTrials(0),
 208      fTrialsRun(0),
 209      fQ(0.),
 210      fX1(0.),
 211      fX2(0.),
 212      fEventTime(0.),
 213      fInteractionRate(0.),
 214      fTimeWindow(0.),
 215      fCurSubEvent(0),
 216      fEventsTime(0),
 217      fNev(0),
 218      fFlavorSelect(0),
 219      fXsection(0.),
 220      fPythia(0),
 221      fPtHardMin(0.),
 222      fPtHardMax(1.e4),
 223      fYHardMin(-1.e10),
 224      fYHardMax(1.e10),
 225      fGinit(kTRUE),
 226      fGfinal(kTRUE),
 227      fHadronisation(kTRUE),
 228      fNpartons(0),
 229      fReadFromFile(kFALSE),
 230      fQuench(kFALSE),
 231      fPtKick(1.),
 232      fFullEvent(kTRUE),
 233      fDecayer(new AliDecayerPythia()),
 234      fDebugEventFirst(-1),
 235      fDebugEventLast(-1),
 236      fEtMinJet(0.),
 237      fEtMaxJet(1.e4),
 238      fEtaMinJet(-20.),
 239      fEtaMaxJet(20.),
 240      fPhiMinJet(0.),
 241      fPhiMaxJet(2.* TMath::Pi()),
 242      fJetReconstruction(kCell),
 243      fEtaMinGamma(-20.),
 244      fEtaMaxGamma(20.),
 245      fPhiMinGamma(0.),
 246      fPhiMaxGamma(2. * TMath::Pi()),
 247      fPycellEtaMax(2.),
 248      fPycellNEta(274),
 249      fPycellNPhi(432),
 250      fPycellThreshold(0.),
 251      fPycellEtSeed(4.),
 252      fPycellMinEtJet(10.),
 253      fPycellMaxRadius(1.),
 254      fStackFillOpt(kFlavorSelection),
 255      fFeedDownOpt(kTRUE),
 256      fFragmentation(kTRUE),
 257      fSetNuclei(kFALSE),
 258      fNewMIS(kFALSE),
 259      fHFoff(kFALSE),
 260      fTriggerParticle(0),
 261      fTriggerEta(0.9),
 262      fCountMode(kCountAll),
 263      fHeader(0),
 264      fRL(0),
 265      fFileName(0)
 266 {
 267 // copy constructor
 268     Pythia.Copy(*this);
 269 }
 270
 271 AliGenPythia::~AliGenPythia()
 272 {
 273 // Destructor
 274   if(fEventsTime) delete fEventsTime;
 275 }
 276
 277 void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
 278 {
 279 // Generate pileup using user specified rate
 280     fInteractionRate = rate;
 281     fTimeWindow = timewindow;
 282     GeneratePileup();
 283 }
 284
 285 void AliGenPythia::GeneratePileup()
 286 {
 287 // Generate sub events time for pileup
 288     fEventsTime = 0;
 289     if(fInteractionRate == 0.) {
 290       Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
 291       return;
 292     }
 293
 294     Int_t npart = NumberParticles();
 295     if(npart < 0) {
 296       Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
 297       return;
 298     }
 299
 300     if(fEventsTime) delete fEventsTime;
 301     fEventsTime = new TArrayF(npart);
 302     TArrayF &array = *fEventsTime;
 303     for(Int_t ipart = 0; ipart < npart; ipart++)
 304       array[ipart] = 0.;
 305
 306     Float_t eventtime = 0.;
 307     while(1)
 308       {
 309         eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
 310         if(eventtime > fTimeWindow) break;
 311         array.Set(array.GetSize()+1);
 312         array[array.GetSize()-1] = eventtime;
 313       }
 314
 315     eventtime = 0.;
 316     while(1)
 317       {
 318         eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
 319         if(TMath::Abs(eventtime) > fTimeWindow) break;
 320         array.Set(array.GetSize()+1);
 321         array[array.GetSize()-1] = eventtime;
 322       }
 323
 324     SetNumberParticles(fEventsTime->GetSize());
 325 }
 326
 327 void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
 328                                        Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
 329 {
 330 // Set pycell parameters
 331     fPycellEtaMax    =  etamax;
 332     fPycellNEta      =  neta;
 333     fPycellNPhi      =  nphi;
 334     fPycellThreshold =  thresh;
 335     fPycellEtSeed    =  etseed;
 336     fPycellMinEtJet  =  minet;
 337     fPycellMaxRadius =  r;
 338 }
 339
 340
 341
 342 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
 343 {
 344   // Set a range of event numbers, for which a table
 345   // of generated particle will be printed
 346   fDebugEventFirst = eventFirst;
 347   fDebugEventLast  = eventLast;
 348   if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
 349 }
 350
 351 void AliGenPythia::Init()
 352 {
 353 // Initialisation
 354
 355     SetMC(AliPythia::Instance());
 356     fPythia=(AliPythia*) fMCEvGen;
 357
 358 //
 359     fParentWeight=1./Float_t(fNpart);
 360 //
 361 //  Forward Paramters to the AliPythia object
 362     fDecayer->SetForceDecay(fForceDecay);
 363     fDecayer->Init();
 364
 365
 366     fPythia->SetCKIN(3,fPtHardMin);
 367     fPythia->SetCKIN(4,fPtHardMax);
 368     fPythia->SetCKIN(7,fYHardMin);
 369     fPythia->SetCKIN(8,fYHardMax);
 370
 371     if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);
 372     // Fragmentation?
 373     if (fFragmentation) {
 374       fPythia->SetMSTP(111,1);
 375     } else {
 376       fPythia->SetMSTP(111,0);
 377     }
 378
 379
 380 //  initial state radiation
 381     fPythia->SetMSTP(61,fGinit);
 382 //  final state radiation
 383     fPythia->SetMSTP(71,fGfinal);
 384 //  pt - kick
 385     if (fPtKick > 0.) {
 386         fPythia->SetMSTP(91,1);
 387         fPythia->SetPARP(91,fPtKick);
 388     } else {
 389         fPythia->SetMSTP(91,0);
 390     }
 391
 392
 393     if (fReadFromFile) {
 394         fRL  =  AliRunLoader::Open(fFileName, "Partons");
 395         fRL->LoadKinematics();
 396         fRL->LoadHeader();
 397     } else {
 398         fRL = 0x0;
 399     }
 400 // Switch off Heavy Flavors on request
 401     if (fHFoff) {
 402         fPythia->SetMSTP(58, 3);
 403         fPythia->SetMSTJ(45, 3);
 404         for (Int_t i = 156; i <= 160; i++) fPythia->SetMDME(i, 1, 0);
 405     }
 406  //
 407     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
 408
 409 //  Parent and Children Selection
 410     switch (fProcess)
 411     {
 412     case kPyOldUEQ2ordered:
 413     case kPyOldUEQ2ordered2:
 414     case kPyOldPopcorn:
 415       break;
 416     case kPyCharm:
 417     case kPyCharmUnforced:
 418     case kPyCharmPbPbMNR:
 419     case kPyCharmpPbMNR:
 420     case kPyCharmppMNR:
 421     case kPyCharmppMNRwmi:
 422         fParentSelect[0] =   411;
 423         fParentSelect[1] =   421;
 424         fParentSelect[2] =   431;
 425         fParentSelect[3] =  4122;
 426         fFlavorSelect    =  4;
 427         break;
 428     case kPyD0PbPbMNR:
 429     case kPyD0pPbMNR:
 430     case kPyD0ppMNR:
 431         fParentSelect[0] =   421;
 432         fFlavorSelect    =   4;
 433         break;
 434     case kPyDPlusPbPbMNR:
 435     case kPyDPluspPbMNR:
 436     case kPyDPlusppMNR:
 437         fParentSelect[0] =   411;
 438         fFlavorSelect    =   4;
 439         break;
 440     case kPyDPlusStrangePbPbMNR:
 441     case kPyDPlusStrangepPbMNR:
 442     case kPyDPlusStrangeppMNR:
 443         fParentSelect[0] =   431;
 444         fFlavorSelect    =   4;
 445         break;
 446     case kPyBeauty:
 447     case kPyBeautyPbPbMNR:
 448     case kPyBeautypPbMNR:
 449     case kPyBeautyppMNR:
 450     case kPyBeautyppMNRwmi:
 451         fParentSelect[0]=  511;
 452         fParentSelect[1]=  521;
 453         fParentSelect[2]=  531;
 454         fParentSelect[3]= 5122;
 455         fParentSelect[4]= 5132;
 456         fParentSelect[5]= 5232;
 457         fParentSelect[6]= 5332;
 458         fFlavorSelect   = 5;
 459         break;
 460     case kPyBeautyUnforced:
 461         fParentSelect[0] =  511;
 462         fParentSelect[1] =  521;
 463         fParentSelect[2] =  531;
 464         fParentSelect[3] = 5122;
 465         fParentSelect[4] = 5132;
 466         fParentSelect[5] = 5232;
 467         fParentSelect[6] = 5332;
 468         fFlavorSelect    = 5;
 469         break;
 470     case kPyJpsiChi:
 471     case kPyJpsi:
 472         fParentSelect[0] = 443;
 473         break;
 474     case kPyMb:
 475     case kPyMbNonDiffr:
 476     case kPyMbMSEL1:
 477     case kPyJets:
 478     case kPyDirectGamma:
 479         break;
 480     case kPyW:
 481     case kPyZ:
 482         break;
 483     }
 484 //
 485 //
 486 //  JetFinder for Trigger
 487 //
 488 //  Configure detector (EMCAL like)
 489 //
 490     fPythia->SetPARU(51, fPycellEtaMax);
 491     fPythia->SetMSTU(51, fPycellNEta);
 492     fPythia->SetMSTU(52, fPycellNPhi);
 493 //
 494 //  Configure Jet Finder
 495 //
 496     fPythia->SetPARU(58,  fPycellThreshold);
 497     fPythia->SetPARU(52,  fPycellEtSeed);
 498     fPythia->SetPARU(53,  fPycellMinEtJet);
 499     fPythia->SetPARU(54,  fPycellMaxRadius);
 500     fPythia->SetMSTU(54,  2);
 501 //
 502 //  This counts the total number of calls to Pyevnt() per run.
 503     fTrialsRun = 0;
 504     fQ         = 0.;
 505     fX1        = 0.;
 506     fX2        = 0.;
 507     fNev       = 0 ;
 508 //
 509 //
 510 //
 511     AliGenMC::Init();
 512 //
 513 //
 514 //
 515     if (fSetNuclei) {
 516         fDyBoost = 0;
 517         Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
 518     }
 519
 520     if (fQuench) {
 521         fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
 522     }
 523 }
 524
 525 void AliGenPythia::Generate()
 526 {
 527 // Generate one event
 528
 529     fDecayer->ForceDecay();
 530
 531     Float_t polar[3]   =   {0,0,0};
 532     Float_t origin[3]  =   {0,0,0};
 533     Float_t p[4];
 534 //  converts from mm/c to s
 535     const Float_t kconv=0.001/2.999792458e8;
 536 //
 537     Int_t nt=0;
 538     Int_t jev=0;
 539     Int_t j, kf;
 540     fTrials=0;
 541     fEventTime = 0.;
 542
 543
 544
 545     //  Set collision vertex position
 546     if (fVertexSmear == kPerEvent) Vertex();
 547
 548 //  event loop
 549     while(1)
 550     {
 551 //
 552 // Produce event
 553 //
 554 //
 555 // Switch hadronisation off
 556 //
 557         fPythia->SetMSTJ(1, 0);
 558 //
 559 // Either produce new event or read partons from file
 560 //
 561         if (!fReadFromFile) {
 562             if (!fNewMIS) {
 563                 fPythia->Pyevnt();
 564             } else {
 565                 fPythia->Pyevnw();
 566             }
 567             fNpartons = fPythia->GetN();
 568         } else {
 569             printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
 570             fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
 571             fPythia->SetN(0);
 572             LoadEvent(fRL->Stack(), 0 , 1);
 573             fPythia->Pyedit(21);
 574         }
 575
 576 //
 577 //  Run quenching routine
 578 //
 579         if (fQuench == 1) {
 580             fPythia->Quench();
 581         } else if (fQuench == 2){
 582             fPythia->Pyquen(208., 0, 0.);
 583         }
 584 //
 585 // Switch hadronisation on
 586 //
 587         fPythia->SetMSTJ(1, 1);
 588 //
 589 // .. and perform hadronisation
 590 //      printf("Calling hadronisation %d\n", fPythia->GetN());
 591         fPythia->Pyexec();
 592         fTrials++;
 593         fPythia->ImportParticles(fParticles,"All");
 594         Boost();
 595 //
 596 //
 597 //
 598         Int_t i;
 599
 600
 601         Int_t np = fParticles->GetEntriesFast();
 602
 603         if (np == 0) continue;
 604 //
 605
 606 //
 607         Int_t* pParent   = new Int_t[np];
 608         Int_t* pSelected = new Int_t[np];
 609         Int_t* trackIt   = new Int_t[np];
 610         for (i = 0; i < np; i++) {
 611             pParent[i]   = -1;
 612             pSelected[i] =  0;
 613             trackIt[i]   =  0;
 614         }
 615
 616         Int_t nc = 0;        // Total n. of selected particles
 617         Int_t nParents = 0;  // Selected parents
 618         Int_t nTkbles = 0;   // Trackable particles
 619         if (fProcess != kPyMb && fProcess != kPyJets &&
 620             fProcess != kPyDirectGamma &&
 621             fProcess != kPyMbNonDiffr  &&
 622             fProcess != kPyMbMSEL1     &&
 623             fProcess != kPyW && fProcess != kPyZ &&
 624             fProcess != kPyCharmppMNRwmi && fProcess != kPyBeautyppMNRwmi) {
 625
 626             for (i = 0; i < np; i++) {
 627                 TParticle* iparticle = (TParticle *) fParticles->At(i);
 628                 Int_t ks = iparticle->GetStatusCode();
 629                 kf = CheckPDGCode(iparticle->GetPdgCode());
 630 // No initial state partons
 631                 if (ks==21) continue;
 632 //
 633 // Heavy Flavor Selection
 634 //
 635                 // quark ?
 636                 kf = TMath::Abs(kf);
 637                 Int_t kfl = kf;
 638                 // Resonance
 639
 640                 if (kfl > 100000) kfl %= 100000;
 641                 if (kfl > 10000)  kfl %= 10000;
 642                 // meson ?
 643                 if  (kfl > 10) kfl/=100;
 644                 // baryon
 645                 if (kfl > 10) kfl/=10;
 646                 Int_t ipa = iparticle->GetFirstMother()-1;
 647                 Int_t kfMo = 0;
 648 //
 649 // Establish mother daughter relation between heavy quarks and mesons
 650 //
 651                 if (kf >= fFlavorSelect && kf <= 6) {
 652                     Int_t idau = iparticle->GetFirstDaughter() - 1;
 653                     if (idau > -1) {
 654                         TParticle* daughter = (TParticle *) fParticles->At(idau);
 655                         Int_t pdgD = daughter->GetPdgCode();
 656                         if (pdgD == 91 || pdgD == 92) {
 657                             Int_t jmin = daughter->GetFirstDaughter() - 1;
 658                             Int_t jmax = daughter->GetLastDaughter()  - 1;
 659                             for (Int_t j = jmin; j <= jmax; j++)
 660                                 ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
 661                         } // is string or cluster
 662                     } // has daughter
 663                 } // heavy quark
 664
 665
 666                 if (ipa > -1) {
 667                     TParticle *  mother = (TParticle *) fParticles->At(ipa);
 668                     kfMo = TMath::Abs(mother->GetPdgCode());
 669                 }
 670
 671                 // What to keep in Stack?
 672                 Bool_t flavorOK = kFALSE;
 673                 Bool_t selectOK = kFALSE;
 674                 if (fFeedDownOpt) {
 675                     if (kfl >= fFlavorSelect) flavorOK = kTRUE;
 676                 } else {
 677                     if (kfl > fFlavorSelect) {
 678                         nc = -1;
 679                         break;
 680                     }
 681                     if (kfl == fFlavorSelect) flavorOK = kTRUE;
 682                 }
 683                 switch (fStackFillOpt) {
 684                 case kFlavorSelection:
 685                     selectOK = kTRUE;
 686                     break;
 687                 case kParentSelection:
 688                     if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
 689                     break;
 690                 }
 691                 if (flavorOK && selectOK) {
 692 //
 693 // Heavy flavor hadron or quark
 694 //
 695 // Kinematic seletion on final state heavy flavor mesons
 696                     if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
 697                     {
 698                         continue;
 699                     }
 700                     pSelected[i] = 1;
 701                     if (ParentSelected(kf)) ++nParents; // Update parent count
 702 //                  printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
 703                 } else {
 704 // Kinematic seletion on decay products
 705                     if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
 706                         && !KinematicSelection(iparticle, 1))
 707                     {
 708                         continue;
 709                     }
 710 //
 711 // Decay products
 712 // Select if mother was selected and is not tracked
 713
 714                     if (pSelected[ipa] &&
 715                         !trackIt[ipa]  &&     // mother will be  tracked ?
 716                         kfMo !=  5 &&         // mother is b-quark, don't store fragments
 717                         kfMo !=  4 &&         // mother is c-quark, don't store fragments
 718                         kf   != 92)           // don't store string
 719                     {
 720 //
 721 // Semi-stable or de-selected: diselect decay products:
 722 //
 723 //
 724                         if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
 725                         {
 726                             Int_t ipF = iparticle->GetFirstDaughter();
 727                             Int_t ipL = iparticle->GetLastDaughter();
 728                             if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
 729                         }
 730 //                      printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
 731                         pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
 732                     }
 733                 }
 734                 if (pSelected[i] == -1) pSelected[i] = 0;
 735                 if (!pSelected[i]) continue;
 736                 // Count quarks only if you did not include fragmentation
 737                 if (fFragmentation && kf <= 10) continue;
 738
 739                 nc++;
 740 // Decision on tracking
 741                 trackIt[i] = 0;
 742 //
 743 // Track final state particle
 744                 if (ks == 1) trackIt[i] = 1;
 745 // Track semi-stable particles
 746                 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
 747 // Track particles selected by process if undecayed.
 748                 if (fForceDecay == kNoDecay) {
 749                     if (ParentSelected(kf)) trackIt[i] = 1;
 750                 } else {
 751                     if (ParentSelected(kf)) trackIt[i] = 0;
 752                 }
 753                 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
 754 //
 755 //
 756
 757             } // particle selection loop
 758             if (nc > 0) {
 759                 for (i = 0; i<np; i++) {
 760                     if (!pSelected[i]) continue;
 761                     TParticle *  iparticle = (TParticle *) fParticles->At(i);
 762                     kf = CheckPDGCode(iparticle->GetPdgCode());
 763                     Int_t ks = iparticle->GetStatusCode();
 764                     p[0] = iparticle->Px();
 765                     p[1] = iparticle->Py();
 766                     p[2] = iparticle->Pz();
 767                     p[3] = iparticle->Energy();
 768
 769                     origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 770                     origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 771                     origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 772
 773                     Float_t tof   = kconv*iparticle->T();
 774                     Int_t ipa     = iparticle->GetFirstMother()-1;
 775                     Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
 776
 777                     PushTrack(fTrackIt*trackIt[i], iparent, kf,
 778                               p[0], p[1], p[2], p[3],
 779                               origin[0], origin[1], origin[2], tof,
 780                               polar[0], polar[1], polar[2],
 781                               kPPrimary, nt, 1., ks);
 782                     pParent[i] = nt;
 783                     KeepTrack(nt);
 784                 } //  PushTrack loop
 785             }
 786         } else {
 787             nc = GenerateMB();
 788         } // mb ?
 789
 790         GetSubEventTime();
 791
 792         delete[] pParent;
 793         delete[] pSelected;
 794         delete[] trackIt;
 795
 796         if (nc > 0) {
 797           switch (fCountMode) {
 798           case kCountAll:
 799             // printf(" Count all \n");
 800             jev += nc;
 801             break;
 802           case kCountParents:
 803             // printf(" Count parents \n");
 804             jev += nParents;
 805             break;
 806           case kCountTrackables:
 807             // printf(" Count trackable \n");
 808             jev += nTkbles;
 809             break;
 810           }
 811             if (jev >= fNpart || fNpart == -1) {
 812                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 813
 814                 fQ  += fPythia->GetVINT(51);
 815                 fX1 += fPythia->GetVINT(41);
 816                 fX2 += fPythia->GetVINT(42);
 817                 fTrialsRun += fTrials;
 818                 fNev++;
 819                 MakeHeader();
 820                 break;
 821             }
 822         }
 823     } // event loop
 824     SetHighWaterMark(nt);
 825 //  adjust weight due to kinematic selection
 826 //    AdjustWeights();
 827 //  get cross-section
 828     fXsection=fPythia->GetPARI(1);
 829 }
 830
 831 Int_t  AliGenPythia::GenerateMB()
 832 {
 833 //
 834 // Min Bias selection and other global selections
 835 //
 836     Int_t i, kf, nt, iparent;
 837     Int_t nc = 0;
 838     Float_t p[4];
 839     Float_t polar[3]   =   {0,0,0};
 840     Float_t origin[3]  =   {0,0,0};
 841 //  converts from mm/c to s
 842     const Float_t kconv=0.001/2.999792458e8;
 843
 844
 845
 846     Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
 847
 848
 849
 850     Int_t* pParent = new Int_t[np];
 851     for (i=0; i< np; i++) pParent[i] = -1;
 852     if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
 853         TParticle* jet1 = (TParticle *) fParticles->At(6);
 854         TParticle* jet2 = (TParticle *) fParticles->At(7);
 855         if (!CheckTrigger(jet1, jet2)) {
 856           delete [] pParent;
 857           return 0;
 858         }
 859     }
 860
 861     if (fTriggerParticle) {
 862         Bool_t triggered = kFALSE;
 863         for (i = 0; i < np; i++) {
 864             TParticle *  iparticle = (TParticle *) fParticles->At(i);
 865             kf = CheckPDGCode(iparticle->GetPdgCode());
 866             if (kf != fTriggerParticle) continue;
 867             if (iparticle->Pt() == 0.) continue;
 868             if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
 869             triggered = kTRUE;
 870             break;
 871         }
 872         if (!triggered) {
 873           delete [] pParent;
 874           return 0;
 875         }
 876     }
 877
 878
 879     // Check if there is a ccbar or bbbar pair with at least one of the two
 880     // in fYMin < y < fYMax
 881     if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
 882       TParticle *hvq;
 883       Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
 884       Float_t yQ;
 885       Int_t   pdgQ;
 886       for(i=0; i<np; i++) {
 887         hvq = (TParticle*)fParticles->At(i);
 888         pdgQ = hvq->GetPdgCode();
 889         if(TMath::Abs(pdgQ) != fFlavorSelect) continue;
 890         if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
 891         yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/
 892                             (hvq->Energy()-hvq->Pz()+1.e-13));
 893         if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE;
 894       }
 895       if (!theQ || !theQbar || !inYcut) {
 896         delete[] pParent;
 897         return 0;
 898       }
 899     }
 900
 901     //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
 902     if ( (fProcess == kPyW || fProcess == kPyZ || fProcess == kPyMb || fProcess == kPyMbNonDiffr)
 903          && (fCutOnChild == 1) ) {
 904       if ( !CheckKinematicsOnChild() ) {
 905         delete[] pParent;
 906         return 0;
 907       }
 908     }
 909
 910
 911     for (i = 0; i < np; i++) {
 912         Int_t trackIt = 0;
 913         TParticle *  iparticle = (TParticle *) fParticles->At(i);
 914         kf = CheckPDGCode(iparticle->GetPdgCode());
 915         Int_t ks = iparticle->GetStatusCode();
 916         Int_t km = iparticle->GetFirstMother();
 917         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
 918             (ks != 1) ||
 919             (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
 920             nc++;
 921             if (ks == 1) trackIt = 1;
 922             Int_t ipa = iparticle->GetFirstMother()-1;
 923
 924             iparent = (ipa > -1) ? pParent[ipa] : -1;
 925
 926 //
 927 // store track information
 928             p[0] = iparticle->Px();
 929             p[1] = iparticle->Py();
 930             p[2] = iparticle->Pz();
 931             p[3] = iparticle->Energy();
 932
 933
 934             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 935             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 936             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 937
 938             Float_t tof = fEventTime + kconv * iparticle->T();
 939
 940             PushTrack(fTrackIt*trackIt, iparent, kf,
 941                       p[0], p[1], p[2], p[3],
 942                       origin[0], origin[1], origin[2], tof,
 943                       polar[0], polar[1], polar[2],
 944                       kPPrimary, nt, 1., ks);
 945             //
 946             // Special Treatment to store color-flow
 947             //
 948             if (ks == 3 || ks == 13 || ks == 14) {
 949                 TParticle* particle = 0;
 950                 if (fStack) {
 951                     particle = fStack->Particle(nt);
 952                 } else {
 953                     particle = gAlice->Stack()->Particle(nt);
 954                 }
 955                 particle->SetFirstDaughter(fPythia->GetK(2, i));
 956                 particle->SetLastDaughter(fPythia->GetK(3, i));
 957             }
 958
 959             KeepTrack(nt);
 960             pParent[i] = nt;
 961             SetHighWaterMark(nt);
 962
 963         } // select particle
 964     } // particle loop
 965
 966     delete[] pParent;
 967
 968     return 1;
 969 }
 970
 971
 972 void AliGenPythia::FinishRun()
 973 {
 974 // Print x-section summary
 975     fPythia->Pystat(1);
 976
 977     if (fNev > 0.) {
 978         fQ  /= fNev;
 979         fX1 /= fNev;
 980         fX2 /= fNev;
 981     }
 982
 983     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
 984     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
 985 }
 986
 987 void AliGenPythia::AdjustWeights() const
 988 {
 989 // Adjust the weights after generation of all events
 990 //
 991     if (gAlice) {
 992         TParticle *part;
 993         Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
 994         for (Int_t i=0; i<ntrack; i++) {
 995             part= gAlice->GetMCApp()->Particle(i);
 996             part->SetWeight(part->GetWeight()*fKineBias);
 997         }
 998     }
 999 }
1000
1001 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
1002 {
1003 // Treat protons as inside nuclei with mass numbers a1 and a2
1004
1005     fAProjectile = a1;
1006     fATarget     = a2;
1007     fSetNuclei   = kTRUE;
1008 }
1009
1010
1011 void AliGenPythia::MakeHeader()
1012 {
1013 //
1014 // Make header for the simulated event
1015 //
1016   if (gAlice) {
1017     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1018         gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1019   }
1020
1021 // Builds the event header, to be called after each event
1022     if (fHeader) delete fHeader;
1023     fHeader = new AliGenPythiaEventHeader("Pythia");
1024 //
1025 // Event type
1026     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1027 //
1028 // Number of trials
1029     ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1030 //
1031 // Event Vertex
1032     fHeader->SetPrimaryVertex(fVertex);
1033 //
1034 // Jets that have triggered
1035
1036     if (fProcess == kPyJets)
1037     {
1038         Int_t ntrig, njet;
1039         Float_t jets[4][10];
1040         GetJets(njet, ntrig, jets);
1041
1042
1043         for (Int_t i = 0; i < ntrig; i++) {
1044             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1045                                                         jets[3][i]);
1046         }
1047     }
1048 //
1049 // Copy relevant information from external header, if present.
1050 //
1051     Float_t uqJet[4];
1052
1053     if (fRL) {
1054         AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1055         for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1056         {
1057             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
1058
1059
1060             exHeader->TriggerJet(i, uqJet);
1061             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1062         }
1063     }
1064 //
1065 // Store quenching parameters
1066 //
1067     if (fQuench){
1068         Double_t z[4];
1069         Double_t xp, yp;
1070         if (fQuench == 1) {
1071             // Pythia::Quench()
1072             fPythia->GetQuenchingParameters(xp, yp, z);
1073         } else {
1074             // Pyquen
1075             Double_t r1 = PARIMP.rb1;
1076             Double_t r2 = PARIMP.rb2;
1077             Double_t b  = PARIMP.b1;
1078             Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1079             Double_t phi = PARIMP.psib1;
1080             xp = r * TMath::Cos(phi);
1081             yp = r * TMath::Sin(phi);
1082
1083         }
1084             ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1085             ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1086         }
1087 //
1088 // Store pt^hard
1089     ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1090 //
1091 //  Pass header
1092 //
1093     AddHeader(fHeader);
1094 }
1095
1096 void AliGenPythia::AddHeader(AliGenEventHeader* header)
1097 {
1098     // Add header to container or runloader
1099     if (fContainer) {
1100         fContainer->AddHeader(header);
1101     } else {
1102         AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
1103     }
1104 }
1105
1106
1107 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
1108 {
1109 // Check the kinematic trigger condition
1110 //
1111     Double_t eta[2];
1112     eta[0] = jet1->Eta();
1113     eta[1] = jet2->Eta();
1114     Double_t phi[2];
1115     phi[0] = jet1->Phi();
1116     phi[1] = jet2->Phi();
1117     Int_t    pdg[2];
1118     pdg[0] = jet1->GetPdgCode();
1119     pdg[1] = jet2->GetPdgCode();
1120     Bool_t   triggered = kFALSE;
1121
1122     if (fProcess == kPyJets) {
1123         Int_t njets = 0;
1124         Int_t ntrig = 0;
1125         Float_t jets[4][10];
1126 //
1127 // Use Pythia clustering on parton level to determine jet axis
1128 //
1129         GetJets(njets, ntrig, jets);
1130
1131         if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1132 //
1133     } else {
1134         Int_t ij = 0;
1135         Int_t ig = 1;
1136         if (pdg[0] == kGamma) {
1137             ij = 1;
1138             ig = 0;
1139         }
1140         //Check eta range first...
1141         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
1142             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1143         {
1144             //Eta is okay, now check phi range
1145             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
1146                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1147             {
1148                 triggered = kTRUE;
1149             }
1150         }
1151     }
1152     return triggered;
1153 }
1154
1155
1156
1157 Bool_t AliGenPythia::CheckKinematicsOnChild(){
1158 //
1159 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1160 //
1161     Bool_t checking = kFALSE;
1162     Int_t j, kcode, ks, km;
1163     Int_t nPartAcc = 0; //number of particles in the acceptance range
1164     Int_t numberOfAcceptedParticles = 1;
1165     if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1166     Int_t npart = fParticles->GetEntriesFast();
1167
1168     for (j = 0; j<npart; j++) {
1169         TParticle *  jparticle = (TParticle *) fParticles->At(j);
1170         kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1171         ks = jparticle->GetStatusCode();
1172         km = jparticle->GetFirstMother();
1173
1174         if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
1175             nPartAcc++;
1176         }
1177         if( numberOfAcceptedParticles <= nPartAcc){
1178           checking = kTRUE;
1179           break;
1180         }
1181     }
1182
1183     return checking;
1184 }
1185
1186
1187 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
1188 {
1189 // Assignment operator
1190     rhs.Copy(*this);
1191     return *this;
1192 }
1193
1194 void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1195 {
1196 //
1197 // Load event into Pythia Common Block
1198 //
1199
1200     Int_t npart = stack -> GetNprimary();
1201     Int_t n0 = 0;
1202
1203     if (!flag) {
1204         (fPythia->GetPyjets())->N = npart;
1205     } else {
1206         n0 = (fPythia->GetPyjets())->N;
1207         (fPythia->GetPyjets())->N = n0 + npart;
1208     }
1209
1210
1211     for (Int_t part = 0; part < npart; part++) {
1212         TParticle *mPart = stack->Particle(part);
1213
1214         Int_t kf     =  mPart->GetPdgCode();
1215         Int_t ks     =  mPart->GetStatusCode();
1216         Int_t idf    =  mPart->GetFirstDaughter();
1217         Int_t idl    =  mPart->GetLastDaughter();
1218
1219         if (reHadr) {
1220             if (ks == 11 || ks == 12) {
1221                 ks  -= 10;
1222                 idf  = -1;
1223                 idl  = -1;
1224             }
1225         }
1226
1227         Float_t px = mPart->Px();
1228         Float_t py = mPart->Py();
1229         Float_t pz = mPart->Pz();
1230         Float_t e  = mPart->Energy();
1231         Float_t m  = mPart->GetCalcMass();
1232
1233
1234         (fPythia->GetPyjets())->P[0][part+n0] = px;
1235         (fPythia->GetPyjets())->P[1][part+n0] = py;
1236         (fPythia->GetPyjets())->P[2][part+n0] = pz;
1237         (fPythia->GetPyjets())->P[3][part+n0] = e;
1238         (fPythia->GetPyjets())->P[4][part+n0] = m;
1239
1240         (fPythia->GetPyjets())->K[1][part+n0] = kf;
1241         (fPythia->GetPyjets())->K[0][part+n0] = ks;
1242         (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1243         (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1244         (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1245     }
1246 }
1247
1248
1249 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1250 {
1251 //
1252 //  Calls the Pythia jet finding algorithm to find jets in the current event
1253 //
1254 //
1255 //
1256 //  Save jets
1257     Int_t n     = fPythia->GetN();
1258
1259 //
1260 //  Run Jet Finder
1261     fPythia->Pycell(njets);
1262     Int_t i;
1263     for (i = 0; i < njets; i++) {
1264         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1265         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1266         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1267         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1268
1269         jets[0][i] = px;
1270         jets[1][i] = py;
1271         jets[2][i] = pz;
1272         jets[3][i] = e;
1273     }
1274 }
1275
1276
1277
1278 void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1279 {
1280 //
1281 //  Calls the Pythia clustering algorithm to find jets in the current event
1282 //
1283     Int_t n     = fPythia->GetN();
1284     nJets       = 0;
1285     nJetsTrig   = 0;
1286     if (fJetReconstruction == kCluster) {
1287 //
1288 //  Configure cluster algorithm
1289 //
1290         fPythia->SetPARU(43, 2.);
1291         fPythia->SetMSTU(41, 1);
1292 //
1293 //  Call cluster algorithm
1294 //
1295         fPythia->Pyclus(nJets);
1296 //
1297 //  Loading jets from common block
1298 //
1299     } else {
1300
1301 //
1302 //  Run Jet Finder
1303         fPythia->Pycell(nJets);
1304     }
1305
1306     Int_t i;
1307     for (i = 0; i < nJets; i++) {
1308         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1309         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1310         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1311         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1312         Float_t pt    = TMath::Sqrt(px * px + py * py);
1313         Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);
1314         Float_t theta = TMath::ATan2(pt,pz);
1315         Float_t et    = e * TMath::Sin(theta);
1316         Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
1317         if (
1318             eta > fEtaMinJet && eta < fEtaMaxJet &&
1319             phi > fPhiMinJet && phi < fPhiMaxJet &&
1320             et  > fEtMinJet  && et  < fEtMaxJet
1321             )
1322         {
1323             jets[0][nJetsTrig] = px;
1324             jets[1][nJetsTrig] = py;
1325             jets[2][nJetsTrig] = pz;
1326             jets[3][nJetsTrig] = e;
1327             nJetsTrig++;
1328 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1329         } else {
1330 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1331         }
1332     }
1333 }
1334
1335 void AliGenPythia::GetSubEventTime()
1336 {
1337   // Calculates time of the next subevent
1338   fEventTime = 0.;
1339   if (fEventsTime) {
1340     TArrayF &array = *fEventsTime;
1341     fEventTime = array[fCurSubEvent++];
1342   }
1343   //  printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1344   return;
1345 }
1346
1347 #ifdef never
1348 void AliGenPythia::Streamer(TBuffer &R__b)
1349 {
1350    // Stream an object of class AliGenPythia.
1351
1352    if (R__b.IsReading()) {
1353       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1354       AliGenerator::Streamer(R__b);
1355       R__b >> (Int_t&)fProcess;
1356       R__b >> (Int_t&)fStrucFunc;
1357       R__b >> (Int_t&)fForceDecay;
1358       R__b >> fEnergyCMS;
1359       R__b >> fKineBias;
1360       R__b >> fTrials;
1361       fParentSelect.Streamer(R__b);
1362       fChildSelect.Streamer(R__b);
1363       R__b >> fXsection;
1364 //      (AliPythia::Instance())->Streamer(R__b);
1365       R__b >> fPtHardMin;
1366       R__b >> fPtHardMax;
1367 //      if (fDecayer) fDecayer->Streamer(R__b);
1368    } else {
1369       R__b.WriteVersion(AliGenPythia::IsA());
1370       AliGenerator::Streamer(R__b);
1371       R__b << (Int_t)fProcess;
1372       R__b << (Int_t)fStrucFunc;
1373       R__b << (Int_t)fForceDecay;
1374       R__b << fEnergyCMS;
1375       R__b << fKineBias;
1376       R__b << fTrials;
1377       fParentSelect.Streamer(R__b);
1378       fChildSelect.Streamer(R__b);
1379       R__b << fXsection;
1380 //      R__b << fPythia;
1381       R__b << fPtHardMin;
1382       R__b << fPtHardMax;
1383       //     fDecayer->Streamer(R__b);
1384    }
1385 }
1386 #endif
1387
1388
1389