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