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