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     fPythia->SetPARJ(200, 0.0);
 524
 525     if (fQuench == 3) {
 526         // Nestor's change of the splittings
 527         fPythia->SetPARJ(200, 0.8);
 528         fPythia->SetMSTJ(41, 1);  // QCD radiation only
 529         fPythia->SetMSTJ(42, 2);  // angular ordering
 530         fPythia->SetMSTJ(44, 2);  // option to run alpha_s
 531         fPythia->SetMSTJ(47, 0);  // No correction back to hard scattering element
 532         fPythia->SetMSTJ(50, 0);  // No coherence in first branching
 533         fPythia->SetPARJ(82, 1.); // Cut off for parton showers
 534     }
 535 }
 536
 537 void AliGenPythia::Generate()
 538 {
 539 // Generate one event
 540
 541     fDecayer->ForceDecay();
 542
 543     Float_t polar[3]   =   {0,0,0};
 544     Float_t origin[3]  =   {0,0,0};
 545     Float_t p[4];
 546 //  converts from mm/c to s
 547     const Float_t kconv=0.001/2.999792458e8;
 548 //
 549     Int_t nt=0;
 550     Int_t jev=0;
 551     Int_t j, kf;
 552     fTrials=0;
 553     fEventTime = 0.;
 554
 555
 556
 557     //  Set collision vertex position
 558     if (fVertexSmear == kPerEvent) Vertex();
 559
 560 //  event loop
 561     while(1)
 562     {
 563 //
 564 // Produce event
 565 //
 566 //
 567 // Switch hadronisation off
 568 //
 569         fPythia->SetMSTJ(1, 0);
 570 //
 571 // Either produce new event or read partons from file
 572 //
 573         if (!fReadFromFile) {
 574             if (!fNewMIS) {
 575                 fPythia->Pyevnt();
 576             } else {
 577                 fPythia->Pyevnw();
 578             }
 579             fNpartons = fPythia->GetN();
 580         } else {
 581             printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
 582             fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
 583             fPythia->SetN(0);
 584             LoadEvent(fRL->Stack(), 0 , 1);
 585             fPythia->Pyedit(21);
 586         }
 587
 588 //
 589 //  Run quenching routine
 590 //
 591         if (fQuench == 1) {
 592             fPythia->Quench();
 593         } else if (fQuench == 2){
 594             fPythia->Pyquen(208., 0, 0.);
 595         } else if (fQuench == 3) {
 596             // Quenching is via multiplicative correction of the splittings
 597         }
 598
 599 //
 600 // Switch hadronisation on
 601 //
 602         fPythia->SetMSTJ(1, 1);
 603 //
 604 // .. and perform hadronisation
 605 //      printf("Calling hadronisation %d\n", fPythia->GetN());
 606         fPythia->Pyexec();
 607         fTrials++;
 608         fPythia->ImportParticles(fParticles,"All");
 609         Boost();
 610 //
 611 //
 612 //
 613         Int_t i;
 614
 615
 616         Int_t np = fParticles->GetEntriesFast();
 617
 618         if (np == 0) continue;
 619 //
 620
 621 //
 622         Int_t* pParent   = new Int_t[np];
 623         Int_t* pSelected = new Int_t[np];
 624         Int_t* trackIt   = new Int_t[np];
 625         for (i = 0; i < np; i++) {
 626             pParent[i]   = -1;
 627             pSelected[i] =  0;
 628             trackIt[i]   =  0;
 629         }
 630
 631         Int_t nc = 0;        // Total n. of selected particles
 632         Int_t nParents = 0;  // Selected parents
 633         Int_t nTkbles = 0;   // Trackable particles
 634         if (fProcess != kPyMb && fProcess != kPyJets &&
 635             fProcess != kPyDirectGamma &&
 636             fProcess != kPyMbNonDiffr  &&
 637             fProcess != kPyMbMSEL1     &&
 638             fProcess != kPyW && fProcess != kPyZ &&
 639             fProcess != kPyCharmppMNRwmi && fProcess != kPyBeautyppMNRwmi) {
 640
 641             for (i = 0; i < np; i++) {
 642                 TParticle* iparticle = (TParticle *) fParticles->At(i);
 643                 Int_t ks = iparticle->GetStatusCode();
 644                 kf = CheckPDGCode(iparticle->GetPdgCode());
 645 // No initial state partons
 646                 if (ks==21) continue;
 647 //
 648 // Heavy Flavor Selection
 649 //
 650                 // quark ?
 651                 kf = TMath::Abs(kf);
 652                 Int_t kfl = kf;
 653                 // Resonance
 654
 655                 if (kfl > 100000) kfl %= 100000;
 656                 if (kfl > 10000)  kfl %= 10000;
 657                 // meson ?
 658                 if  (kfl > 10) kfl/=100;
 659                 // baryon
 660                 if (kfl > 10) kfl/=10;
 661                 Int_t ipa = iparticle->GetFirstMother()-1;
 662                 Int_t kfMo = 0;
 663 //
 664 // Establish mother daughter relation between heavy quarks and mesons
 665 //
 666                 if (kf >= fFlavorSelect && kf <= 6) {
 667                     Int_t idau = iparticle->GetFirstDaughter() - 1;
 668                     if (idau > -1) {
 669                         TParticle* daughter = (TParticle *) fParticles->At(idau);
 670                         Int_t pdgD = daughter->GetPdgCode();
 671                         if (pdgD == 91 || pdgD == 92) {
 672                             Int_t jmin = daughter->GetFirstDaughter() - 1;
 673                             Int_t jmax = daughter->GetLastDaughter()  - 1;
 674                             for (Int_t j = jmin; j <= jmax; j++)
 675                                 ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
 676                         } // is string or cluster
 677                     } // has daughter
 678                 } // heavy quark
 679
 680
 681                 if (ipa > -1) {
 682                     TParticle *  mother = (TParticle *) fParticles->At(ipa);
 683                     kfMo = TMath::Abs(mother->GetPdgCode());
 684                 }
 685
 686                 // What to keep in Stack?
 687                 Bool_t flavorOK = kFALSE;
 688                 Bool_t selectOK = kFALSE;
 689                 if (fFeedDownOpt) {
 690                     if (kfl >= fFlavorSelect) flavorOK = kTRUE;
 691                 } else {
 692                     if (kfl > fFlavorSelect) {
 693                         nc = -1;
 694                         break;
 695                     }
 696                     if (kfl == fFlavorSelect) flavorOK = kTRUE;
 697                 }
 698                 switch (fStackFillOpt) {
 699                 case kFlavorSelection:
 700                     selectOK = kTRUE;
 701                     break;
 702                 case kParentSelection:
 703                     if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
 704                     break;
 705                 }
 706                 if (flavorOK && selectOK) {
 707 //
 708 // Heavy flavor hadron or quark
 709 //
 710 // Kinematic seletion on final state heavy flavor mesons
 711                     if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
 712                     {
 713                         continue;
 714                     }
 715                     pSelected[i] = 1;
 716                     if (ParentSelected(kf)) ++nParents; // Update parent count
 717 //                  printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
 718                 } else {
 719 // Kinematic seletion on decay products
 720                     if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
 721                         && !KinematicSelection(iparticle, 1))
 722                     {
 723                         continue;
 724                     }
 725 //
 726 // Decay products
 727 // Select if mother was selected and is not tracked
 728
 729                     if (pSelected[ipa] &&
 730                         !trackIt[ipa]  &&     // mother will be  tracked ?
 731                         kfMo !=  5 &&         // mother is b-quark, don't store fragments
 732                         kfMo !=  4 &&         // mother is c-quark, don't store fragments
 733                         kf   != 92)           // don't store string
 734                     {
 735 //
 736 // Semi-stable or de-selected: diselect decay products:
 737 //
 738 //
 739                         if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
 740                         {
 741                             Int_t ipF = iparticle->GetFirstDaughter();
 742                             Int_t ipL = iparticle->GetLastDaughter();
 743                             if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
 744                         }
 745 //                      printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
 746                         pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
 747                     }
 748                 }
 749                 if (pSelected[i] == -1) pSelected[i] = 0;
 750                 if (!pSelected[i]) continue;
 751                 // Count quarks only if you did not include fragmentation
 752                 if (fFragmentation && kf <= 10) continue;
 753
 754                 nc++;
 755 // Decision on tracking
 756                 trackIt[i] = 0;
 757 //
 758 // Track final state particle
 759                 if (ks == 1) trackIt[i] = 1;
 760 // Track semi-stable particles
 761                 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
 762 // Track particles selected by process if undecayed.
 763                 if (fForceDecay == kNoDecay) {
 764                     if (ParentSelected(kf)) trackIt[i] = 1;
 765                 } else {
 766                     if (ParentSelected(kf)) trackIt[i] = 0;
 767                 }
 768                 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
 769 //
 770 //
 771
 772             } // particle selection loop
 773             if (nc > 0) {
 774                 for (i = 0; i<np; i++) {
 775                     if (!pSelected[i]) continue;
 776                     TParticle *  iparticle = (TParticle *) fParticles->At(i);
 777                     kf = CheckPDGCode(iparticle->GetPdgCode());
 778                     Int_t ks = iparticle->GetStatusCode();
 779                     p[0] = iparticle->Px();
 780                     p[1] = iparticle->Py();
 781                     p[2] = iparticle->Pz();
 782                     p[3] = iparticle->Energy();
 783
 784                     origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 785                     origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 786                     origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 787
 788                     Float_t tof   = kconv*iparticle->T();
 789                     Int_t ipa     = iparticle->GetFirstMother()-1;
 790                     Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
 791
 792                     PushTrack(fTrackIt*trackIt[i], iparent, kf,
 793                               p[0], p[1], p[2], p[3],
 794                               origin[0], origin[1], origin[2], tof,
 795                               polar[0], polar[1], polar[2],
 796                               kPPrimary, nt, 1., ks);
 797                     pParent[i] = nt;
 798                     KeepTrack(nt);
 799                 } //  PushTrack loop
 800             }
 801         } else {
 802             nc = GenerateMB();
 803         } // mb ?
 804
 805         GetSubEventTime();
 806
 807         delete[] pParent;
 808         delete[] pSelected;
 809         delete[] trackIt;
 810
 811         if (nc > 0) {
 812           switch (fCountMode) {
 813           case kCountAll:
 814             // printf(" Count all \n");
 815             jev += nc;
 816             break;
 817           case kCountParents:
 818             // printf(" Count parents \n");
 819             jev += nParents;
 820             break;
 821           case kCountTrackables:
 822             // printf(" Count trackable \n");
 823             jev += nTkbles;
 824             break;
 825           }
 826             if (jev >= fNpart || fNpart == -1) {
 827                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 828
 829                 fQ  += fPythia->GetVINT(51);
 830                 fX1 += fPythia->GetVINT(41);
 831                 fX2 += fPythia->GetVINT(42);
 832                 fTrialsRun += fTrials;
 833                 fNev++;
 834                 MakeHeader();
 835                 break;
 836             }
 837         }
 838     } // event loop
 839     SetHighWaterMark(nt);
 840 //  adjust weight due to kinematic selection
 841 //    AdjustWeights();
 842 //  get cross-section
 843     fXsection=fPythia->GetPARI(1);
 844 }
 845
 846 Int_t  AliGenPythia::GenerateMB()
 847 {
 848 //
 849 // Min Bias selection and other global selections
 850 //
 851     Int_t i, kf, nt, iparent;
 852     Int_t nc = 0;
 853     Float_t p[4];
 854     Float_t polar[3]   =   {0,0,0};
 855     Float_t origin[3]  =   {0,0,0};
 856 //  converts from mm/c to s
 857     const Float_t kconv=0.001/2.999792458e8;
 858
 859
 860
 861     Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
 862
 863
 864
 865     Int_t* pParent = new Int_t[np];
 866     for (i=0; i< np; i++) pParent[i] = -1;
 867     if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
 868         TParticle* jet1 = (TParticle *) fParticles->At(6);
 869         TParticle* jet2 = (TParticle *) fParticles->At(7);
 870         if (!CheckTrigger(jet1, jet2)) {
 871           delete [] pParent;
 872           return 0;
 873         }
 874     }
 875
 876     if (fTriggerParticle) {
 877         Bool_t triggered = kFALSE;
 878         for (i = 0; i < np; i++) {
 879             TParticle *  iparticle = (TParticle *) fParticles->At(i);
 880             kf = CheckPDGCode(iparticle->GetPdgCode());
 881             if (kf != fTriggerParticle) continue;
 882             if (iparticle->Pt() == 0.) continue;
 883             if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
 884             triggered = kTRUE;
 885             break;
 886         }
 887         if (!triggered) {
 888           delete [] pParent;
 889           return 0;
 890         }
 891     }
 892
 893
 894     // Check if there is a ccbar or bbbar pair with at least one of the two
 895     // in fYMin < y < fYMax
 896     if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
 897       TParticle *hvq;
 898       Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
 899       Float_t yQ;
 900       Int_t   pdgQ;
 901       for(i=0; i<np; i++) {
 902         hvq = (TParticle*)fParticles->At(i);
 903         pdgQ = hvq->GetPdgCode();
 904         if(TMath::Abs(pdgQ) != fFlavorSelect) continue;
 905         if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
 906         yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/
 907                             (hvq->Energy()-hvq->Pz()+1.e-13));
 908         if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE;
 909       }
 910       if (!theQ || !theQbar || !inYcut) {
 911         delete[] pParent;
 912         return 0;
 913       }
 914     }
 915
 916     //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
 917     if ( (fProcess == kPyW || fProcess == kPyZ || fProcess == kPyMb || fProcess == kPyMbNonDiffr)
 918          && (fCutOnChild == 1) ) {
 919       if ( !CheckKinematicsOnChild() ) {
 920         delete[] pParent;
 921         return 0;
 922       }
 923     }
 924
 925
 926     for (i = 0; i < np; i++) {
 927         Int_t trackIt = 0;
 928         TParticle *  iparticle = (TParticle *) fParticles->At(i);
 929         kf = CheckPDGCode(iparticle->GetPdgCode());
 930         Int_t ks = iparticle->GetStatusCode();
 931         Int_t km = iparticle->GetFirstMother();
 932         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
 933             (ks != 1) ||
 934             (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
 935             nc++;
 936             if (ks == 1) trackIt = 1;
 937             Int_t ipa = iparticle->GetFirstMother()-1;
 938
 939             iparent = (ipa > -1) ? pParent[ipa] : -1;
 940
 941 //
 942 // store track information
 943             p[0] = iparticle->Px();
 944             p[1] = iparticle->Py();
 945             p[2] = iparticle->Pz();
 946             p[3] = iparticle->Energy();
 947
 948
 949             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 950             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 951             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 952
 953             Float_t tof = fEventTime + kconv * iparticle->T();
 954
 955             PushTrack(fTrackIt*trackIt, iparent, kf,
 956                       p[0], p[1], p[2], p[3],
 957                       origin[0], origin[1], origin[2], tof,
 958                       polar[0], polar[1], polar[2],
 959                       kPPrimary, nt, 1., ks);
 960             //
 961             // Special Treatment to store color-flow
 962             //
 963             if (ks == 3 || ks == 13 || ks == 14) {
 964                 TParticle* particle = 0;
 965                 if (fStack) {
 966                     particle = fStack->Particle(nt);
 967                 } else {
 968                     particle = gAlice->Stack()->Particle(nt);
 969                 }
 970                 particle->SetFirstDaughter(fPythia->GetK(2, i));
 971                 particle->SetLastDaughter(fPythia->GetK(3, i));
 972             }
 973
 974             KeepTrack(nt);
 975             pParent[i] = nt;
 976             SetHighWaterMark(nt);
 977
 978         } // select particle
 979     } // particle loop
 980
 981     delete[] pParent;
 982
 983     return 1;
 984 }
 985
 986
 987 void AliGenPythia::FinishRun()
 988 {
 989 // Print x-section summary
 990     fPythia->Pystat(1);
 991
 992     if (fNev > 0.) {
 993         fQ  /= fNev;
 994         fX1 /= fNev;
 995         fX2 /= fNev;
 996     }
 997
 998     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
 999     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1000 }
1001
1002 void AliGenPythia::AdjustWeights() const
1003 {
1004 // Adjust the weights after generation of all events
1005 //
1006     if (gAlice) {
1007         TParticle *part;
1008         Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1009         for (Int_t i=0; i<ntrack; i++) {
1010             part= gAlice->GetMCApp()->Particle(i);
1011             part->SetWeight(part->GetWeight()*fKineBias);
1012         }
1013     }
1014 }
1015
1016 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
1017 {
1018 // Treat protons as inside nuclei with mass numbers a1 and a2
1019
1020     fAProjectile = a1;
1021     fATarget     = a2;
1022     fSetNuclei   = kTRUE;
1023 }
1024
1025
1026 void AliGenPythia::MakeHeader()
1027 {
1028 //
1029 // Make header for the simulated event
1030 //
1031   if (gAlice) {
1032     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1033         gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1034   }
1035
1036 // Builds the event header, to be called after each event
1037     if (fHeader) delete fHeader;
1038     fHeader = new AliGenPythiaEventHeader("Pythia");
1039 //
1040 // Event type
1041     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1042 //
1043 // Number of trials
1044     ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1045 //
1046 // Event Vertex
1047     fHeader->SetPrimaryVertex(fVertex);
1048 //
1049 // Jets that have triggered
1050
1051     if (fProcess == kPyJets)
1052     {
1053         Int_t ntrig, njet;
1054         Float_t jets[4][10];
1055         GetJets(njet, ntrig, jets);
1056
1057
1058         for (Int_t i = 0; i < ntrig; i++) {
1059             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1060                                                         jets[3][i]);
1061         }
1062     }
1063 //
1064 // Copy relevant information from external header, if present.
1065 //
1066     Float_t uqJet[4];
1067
1068     if (fRL) {
1069         AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1070         for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1071         {
1072             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
1073
1074
1075             exHeader->TriggerJet(i, uqJet);
1076             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1077         }
1078     }
1079 //
1080 // Store quenching parameters
1081 //
1082     if (fQuench){
1083         Double_t z[4];
1084         Double_t xp, yp;
1085         if (fQuench == 1) {
1086             // Pythia::Quench()
1087             fPythia->GetQuenchingParameters(xp, yp, z);
1088         } else {
1089             // Pyquen
1090             Double_t r1 = PARIMP.rb1;
1091             Double_t r2 = PARIMP.rb2;
1092             Double_t b  = PARIMP.b1;
1093             Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1094             Double_t phi = PARIMP.psib1;
1095             xp = r * TMath::Cos(phi);
1096             yp = r * TMath::Sin(phi);
1097
1098         }
1099             ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1100             ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1101         }
1102 //
1103 // Store pt^hard
1104     ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1105 //
1106 //  Pass header
1107 //
1108     AddHeader(fHeader);
1109 }
1110
1111 void AliGenPythia::AddHeader(AliGenEventHeader* header)
1112 {
1113     // Add header to container or runloader
1114     if (fContainer) {
1115         fContainer->AddHeader(header);
1116     } else {
1117         AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
1118     }
1119 }
1120
1121
1122 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
1123 {
1124 // Check the kinematic trigger condition
1125 //
1126     Double_t eta[2];
1127     eta[0] = jet1->Eta();
1128     eta[1] = jet2->Eta();
1129     Double_t phi[2];
1130     phi[0] = jet1->Phi();
1131     phi[1] = jet2->Phi();
1132     Int_t    pdg[2];
1133     pdg[0] = jet1->GetPdgCode();
1134     pdg[1] = jet2->GetPdgCode();
1135     Bool_t   triggered = kFALSE;
1136
1137     if (fProcess == kPyJets) {
1138         Int_t njets = 0;
1139         Int_t ntrig = 0;
1140         Float_t jets[4][10];
1141 //
1142 // Use Pythia clustering on parton level to determine jet axis
1143 //
1144         GetJets(njets, ntrig, jets);
1145
1146         if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1147 //
1148     } else {
1149         Int_t ij = 0;
1150         Int_t ig = 1;
1151         if (pdg[0] == kGamma) {
1152             ij = 1;
1153             ig = 0;
1154         }
1155         //Check eta range first...
1156         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
1157             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1158         {
1159             //Eta is okay, now check phi range
1160             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
1161                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1162             {
1163                 triggered = kTRUE;
1164             }
1165         }
1166     }
1167     return triggered;
1168 }
1169
1170
1171
1172 Bool_t AliGenPythia::CheckKinematicsOnChild(){
1173 //
1174 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1175 //
1176     Bool_t checking = kFALSE;
1177     Int_t j, kcode, ks, km;
1178     Int_t nPartAcc = 0; //number of particles in the acceptance range
1179     Int_t numberOfAcceptedParticles = 1;
1180     if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1181     Int_t npart = fParticles->GetEntriesFast();
1182
1183     for (j = 0; j<npart; j++) {
1184         TParticle *  jparticle = (TParticle *) fParticles->At(j);
1185         kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1186         ks = jparticle->GetStatusCode();
1187         km = jparticle->GetFirstMother();
1188
1189         if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
1190             nPartAcc++;
1191         }
1192         if( numberOfAcceptedParticles <= nPartAcc){
1193           checking = kTRUE;
1194           break;
1195         }
1196     }
1197
1198     return checking;
1199 }
1200
1201
1202 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
1203 {
1204 // Assignment operator
1205     rhs.Copy(*this);
1206     return *this;
1207 }
1208
1209 void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1210 {
1211 //
1212 // Load event into Pythia Common Block
1213 //
1214
1215     Int_t npart = stack -> GetNprimary();
1216     Int_t n0 = 0;
1217
1218     if (!flag) {
1219         (fPythia->GetPyjets())->N = npart;
1220     } else {
1221         n0 = (fPythia->GetPyjets())->N;
1222         (fPythia->GetPyjets())->N = n0 + npart;
1223     }
1224
1225
1226     for (Int_t part = 0; part < npart; part++) {
1227         TParticle *mPart = stack->Particle(part);
1228
1229         Int_t kf     =  mPart->GetPdgCode();
1230         Int_t ks     =  mPart->GetStatusCode();
1231         Int_t idf    =  mPart->GetFirstDaughter();
1232         Int_t idl    =  mPart->GetLastDaughter();
1233
1234         if (reHadr) {
1235             if (ks == 11 || ks == 12) {
1236                 ks  -= 10;
1237                 idf  = -1;
1238                 idl  = -1;
1239             }
1240         }
1241
1242         Float_t px = mPart->Px();
1243         Float_t py = mPart->Py();
1244         Float_t pz = mPart->Pz();
1245         Float_t e  = mPart->Energy();
1246         Float_t m  = mPart->GetCalcMass();
1247
1248
1249         (fPythia->GetPyjets())->P[0][part+n0] = px;
1250         (fPythia->GetPyjets())->P[1][part+n0] = py;
1251         (fPythia->GetPyjets())->P[2][part+n0] = pz;
1252         (fPythia->GetPyjets())->P[3][part+n0] = e;
1253         (fPythia->GetPyjets())->P[4][part+n0] = m;
1254
1255         (fPythia->GetPyjets())->K[1][part+n0] = kf;
1256         (fPythia->GetPyjets())->K[0][part+n0] = ks;
1257         (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1258         (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1259         (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1260     }
1261 }
1262
1263
1264 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1265 {
1266 //
1267 //  Calls the Pythia jet finding algorithm to find jets in the current event
1268 //
1269 //
1270 //
1271 //  Save jets
1272     Int_t n     = fPythia->GetN();
1273
1274 //
1275 //  Run Jet Finder
1276     fPythia->Pycell(njets);
1277     Int_t i;
1278     for (i = 0; i < njets; i++) {
1279         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1280         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1281         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1282         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1283
1284         jets[0][i] = px;
1285         jets[1][i] = py;
1286         jets[2][i] = pz;
1287         jets[3][i] = e;
1288     }
1289 }
1290
1291
1292
1293 void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1294 {
1295 //
1296 //  Calls the Pythia clustering algorithm to find jets in the current event
1297 //
1298     Int_t n     = fPythia->GetN();
1299     nJets       = 0;
1300     nJetsTrig   = 0;
1301     if (fJetReconstruction == kCluster) {
1302 //
1303 //  Configure cluster algorithm
1304 //
1305         fPythia->SetPARU(43, 2.);
1306         fPythia->SetMSTU(41, 1);
1307 //
1308 //  Call cluster algorithm
1309 //
1310         fPythia->Pyclus(nJets);
1311 //
1312 //  Loading jets from common block
1313 //
1314     } else {
1315
1316 //
1317 //  Run Jet Finder
1318         fPythia->Pycell(nJets);
1319     }
1320
1321     Int_t i;
1322     for (i = 0; i < nJets; i++) {
1323         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1324         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1325         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1326         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1327         Float_t pt    = TMath::Sqrt(px * px + py * py);
1328         Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);
1329         Float_t theta = TMath::ATan2(pt,pz);
1330         Float_t et    = e * TMath::Sin(theta);
1331         Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
1332         if (
1333             eta > fEtaMinJet && eta < fEtaMaxJet &&
1334             phi > fPhiMinJet && phi < fPhiMaxJet &&
1335             et  > fEtMinJet  && et  < fEtMaxJet
1336             )
1337         {
1338             jets[0][nJetsTrig] = px;
1339             jets[1][nJetsTrig] = py;
1340             jets[2][nJetsTrig] = pz;
1341             jets[3][nJetsTrig] = e;
1342             nJetsTrig++;
1343 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1344         } else {
1345 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1346         }
1347     }
1348 }
1349
1350 void AliGenPythia::GetSubEventTime()
1351 {
1352   // Calculates time of the next subevent
1353   fEventTime = 0.;
1354   if (fEventsTime) {
1355     TArrayF &array = *fEventsTime;
1356     fEventTime = array[fCurSubEvent++];
1357   }
1358   //  printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1359   return;
1360 }
1361
1362 #ifdef never
1363 void AliGenPythia::Streamer(TBuffer &R__b)
1364 {
1365    // Stream an object of class AliGenPythia.
1366
1367    if (R__b.IsReading()) {
1368       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1369       AliGenerator::Streamer(R__b);
1370       R__b >> (Int_t&)fProcess;
1371       R__b >> (Int_t&)fStrucFunc;
1372       R__b >> (Int_t&)fForceDecay;
1373       R__b >> fEnergyCMS;
1374       R__b >> fKineBias;
1375       R__b >> fTrials;
1376       fParentSelect.Streamer(R__b);
1377       fChildSelect.Streamer(R__b);
1378       R__b >> fXsection;
1379 //      (AliPythia::Instance())->Streamer(R__b);
1380       R__b >> fPtHardMin;
1381       R__b >> fPtHardMax;
1382 //      if (fDecayer) fDecayer->Streamer(R__b);
1383    } else {
1384       R__b.WriteVersion(AliGenPythia::IsA());
1385       AliGenerator::Streamer(R__b);
1386       R__b << (Int_t)fProcess;
1387       R__b << (Int_t)fStrucFunc;
1388       R__b << (Int_t)fForceDecay;
1389       R__b << fEnergyCMS;
1390       R__b << fKineBias;
1391       R__b << fTrials;
1392       fParentSelect.Streamer(R__b);
1393       fChildSelect.Streamer(R__b);
1394       R__b << fXsection;
1395 //      R__b << fPythia;
1396       R__b << fPtHardMin;
1397       R__b << fPtHardMax;
1398       //     fDecayer->Streamer(R__b);
1399    }
1400 }
1401 #endif
1402
1403
1404