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