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