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