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