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 kPyMb:
 465     case kPyMbWithDirectPhoton:
 466     case kPyMbNonDiffr:
 467     case kPyMbMSEL1:
 468     case kPyJets:
 469     case kPyDirectGamma:
 470     case kPyLhwgMb:
 471         break;
 472     case kPyW:
 473     case kPyZ:
 474         break;
 475     }
 476 //
 477 //
 478 //  JetFinder for Trigger
 479 //
 480 //  Configure detector (EMCAL like)
 481 //
 482     fPythia->SetPARU(51, fPycellEtaMax);
 483     fPythia->SetMSTU(51, fPycellNEta);
 484     fPythia->SetMSTU(52, fPycellNPhi);
 485 //
 486 //  Configure Jet Finder
 487 //
 488     fPythia->SetPARU(58,  fPycellThreshold);
 489     fPythia->SetPARU(52,  fPycellEtSeed);
 490     fPythia->SetPARU(53,  fPycellMinEtJet);
 491     fPythia->SetPARU(54,  fPycellMaxRadius);
 492     fPythia->SetMSTU(54,  2);
 493 //
 494 //  This counts the total number of calls to Pyevnt() per run.
 495     fTrialsRun = 0;
 496     fQ         = 0.;
 497     fX1        = 0.;
 498     fX2        = 0.;
 499     fNev       = 0 ;
 500 //
 501 //
 502 //
 503     AliGenMC::Init();
 504 //
 505 //
 506 //
 507     if (fSetNuclei) {
 508         fDyBoost = 0;
 509         Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
 510     }
 511
 512     fPythia->SetPARJ(200, 0.0);
 513     fPythia->SetPARJ(199, 0.0);
 514     fPythia->SetPARJ(198, 0.0);
 515     fPythia->SetPARJ(197, 0.0);
 516
 517     if (fQuench == 1) {
 518         fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
 519     }
 520
 521     if (fQuench == 3) {
 522         // Nestor's change of the splittings
 523         fPythia->SetPARJ(200, 0.8);
 524         fPythia->SetMSTJ(41, 1);  // QCD radiation only
 525         fPythia->SetMSTJ(42, 2);  // angular ordering
 526         fPythia->SetMSTJ(44, 2);  // option to run alpha_s
 527         fPythia->SetMSTJ(47, 0);  // No correction back to hard scattering element
 528         fPythia->SetMSTJ(50, 0);  // No coherence in first branching
 529         fPythia->SetPARJ(82, 1.); // Cut off for parton showers
 530     } else if (fQuench == 4) {
 531         // Armesto-Cunqueiro-Salgado change of the splittings.
 532         AliFastGlauber* glauber = AliFastGlauber::Instance();
 533         glauber->Init(2);
 534         //read and store transverse almonds corresponding to differnt
 535         //impact parameters.
 536         glauber->SetCentralityClass(0.,0.1);
 537         fPythia->SetPARJ(200, 1.);
 538         fPythia->SetPARJ(198, fQhat);
 539         fPythia->SetPARJ(199, fLength);
 540         fPythia->SetMSTJ(42, 2);  // angular ordering
 541         fPythia->SetMSTJ(44, 2);  // option to run alpha_s
 542         fPythia->SetPARJ(82, 1.); // Cut off for parton showers
 543     }
 544 }
 545
 546 void AliGenPythia::Generate()
 547 {
 548 // Generate one event
 549     if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
 550     fDecayer->ForceDecay();
 551
 552     Float_t polar[3]   =   {0,0,0};
 553     Float_t origin[3]  =   {0,0,0};
 554     Float_t p[4];
 555 //  converts from mm/c to s
 556     const Float_t kconv=0.001/2.999792458e8;
 557 //
 558     Int_t nt=0;
 559     Int_t jev=0;
 560     Int_t j, kf;
 561     fTrials=0;
 562     fEventTime = 0.;
 563
 564
 565
 566     //  Set collision vertex position
 567     if (fVertexSmear == kPerEvent) Vertex();
 568
 569 //  event loop
 570     while(1)
 571     {
 572 //
 573 // Produce event
 574 //
 575 //
 576 // Switch hadronisation off
 577 //
 578         fPythia->SetMSTJ(1, 0);
 579
 580         if (fQuench ==4){
 581             Double_t bimp;
 582             // Quenching comes through medium-modified splitting functions.
 583             AliFastGlauber::Instance()->GetRandomBHard(bimp);
 584             fPythia->SetPARJ(197, bimp);
 585             fImpact = bimp;
 586         }
 587 //
 588 // Either produce new event or read partons from file
 589 //
 590         if (!fReadFromFile) {
 591             if (!fNewMIS) {
 592                 fPythia->Pyevnt();
 593             } else {
 594                 fPythia->Pyevnw();
 595             }
 596             fNpartons = fPythia->GetN();
 597         } else {
 598             printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
 599             fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
 600             fPythia->SetN(0);
 601             LoadEvent(fRL->Stack(), 0 , 1);
 602             fPythia->Pyedit(21);
 603         }
 604
 605 //
 606 //  Run quenching routine
 607 //
 608         if (fQuench == 1) {
 609             fPythia->Quench();
 610         } else if (fQuench == 2){
 611             fPythia->Pyquen(208., 0, 0.);
 612         } else if (fQuench == 3) {
 613             // Quenching is via multiplicative correction of the splittings
 614         }
 615
 616 //
 617 // Switch hadronisation on
 618 //
 619         if (fHadronisation) {
 620             fPythia->SetMSTJ(1, 1);
 621 //
 622 // .. and perform hadronisation
 623 //      printf("Calling hadronisation %d\n", fPythia->GetN());
 624             fPythia->Pyexec();
 625         }
 626         fTrials++;
 627         fPythia->ImportParticles(&fParticles,"All");
 628         if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
 629 //
 630 //
 631 //
 632         Int_t i;
 633
 634         fNprimaries = 0;
 635         Int_t np = fParticles.GetEntriesFast();
 636
 637         if (np == 0) continue;
 638 //
 639
 640 //
 641         Int_t* pParent   = new Int_t[np];
 642         Int_t* pSelected = new Int_t[np];
 643         Int_t* trackIt   = new Int_t[np];
 644         for (i = 0; i < np; i++) {
 645             pParent[i]   = -1;
 646             pSelected[i] =  0;
 647             trackIt[i]   =  0;
 648         }
 649
 650         Int_t nc = 0;        // Total n. of selected particles
 651         Int_t nParents = 0;  // Selected parents
 652         Int_t nTkbles = 0;   // Trackable particles
 653         if (fProcess != kPyMbDefault &&
 654             fProcess != kPyMb &&
 655             fProcess != kPyMbWithDirectPhoton &&
 656             fProcess != kPyJets &&
 657             fProcess != kPyDirectGamma &&
 658             fProcess != kPyMbNonDiffr  &&
 659             fProcess != kPyMbMSEL1     &&
 660             fProcess != kPyW &&
 661             fProcess != kPyZ &&
 662             fProcess != kPyCharmppMNRwmi &&
 663             fProcess != kPyBeautyppMNRwmi &&
 664             fProcess != kPyBeautyJets) {
 665
 666             for (i = 0; i < np; i++) {
 667                 TParticle* iparticle = (TParticle *) fParticles.At(i);
 668                 Int_t ks = iparticle->GetStatusCode();
 669                 kf = CheckPDGCode(iparticle->GetPdgCode());
 670 // No initial state partons
 671                 if (ks==21) continue;
 672 //
 673 // Heavy Flavor Selection
 674 //
 675                 // quark ?
 676                 kf = TMath::Abs(kf);
 677                 Int_t kfl = kf;
 678                 // Resonance
 679
 680                 if (kfl > 100000) kfl %= 100000;
 681                 if (kfl > 10000)  kfl %= 10000;
 682                 // meson ?
 683                 if  (kfl > 10) kfl/=100;
 684                 // baryon
 685                 if (kfl > 10) kfl/=10;
 686                 Int_t ipa = iparticle->GetFirstMother()-1;
 687                 Int_t kfMo = 0;
 688 //
 689 // Establish mother daughter relation between heavy quarks and mesons
 690 //
 691                 if (kf >= fFlavorSelect && kf <= 6) {
 692                     Int_t idau = iparticle->GetFirstDaughter() - 1;
 693                     if (idau > -1) {
 694                         TParticle* daughter = (TParticle *) fParticles.At(idau);
 695                         Int_t pdgD = daughter->GetPdgCode();
 696                         if (pdgD == 91 || pdgD == 92) {
 697                             Int_t jmin = daughter->GetFirstDaughter() - 1;
 698                             Int_t jmax = daughter->GetLastDaughter()  - 1;
 699                             for (Int_t jp = jmin; jp <= jmax; jp++)
 700                                 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
 701                         } // is string or cluster
 702                     } // has daughter
 703                 } // heavy quark
 704
 705
 706                 if (ipa > -1) {
 707                     TParticle *  mother = (TParticle *) fParticles.At(ipa);
 708                     kfMo = TMath::Abs(mother->GetPdgCode());
 709                 }
 710
 711                 // What to keep in Stack?
 712                 Bool_t flavorOK = kFALSE;
 713                 Bool_t selectOK = kFALSE;
 714                 if (fFeedDownOpt) {
 715                     if (kfl >= fFlavorSelect) flavorOK = kTRUE;
 716                 } else {
 717                     if (kfl > fFlavorSelect) {
 718                         nc = -1;
 719                         break;
 720                     }
 721                     if (kfl == fFlavorSelect) flavorOK = kTRUE;
 722                 }
 723                 switch (fStackFillOpt) {
 724                 case kFlavorSelection:
 725                     selectOK = kTRUE;
 726                     break;
 727                 case kParentSelection:
 728                     if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
 729                     break;
 730                 }
 731                 if (flavorOK && selectOK) {
 732 //
 733 // Heavy flavor hadron or quark
 734 //
 735 // Kinematic seletion on final state heavy flavor mesons
 736                     if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
 737                     {
 738                         continue;
 739                     }
 740                     pSelected[i] = 1;
 741                     if (ParentSelected(kf)) ++nParents; // Update parent count
 742 //                  printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
 743                 } else {
 744 // Kinematic seletion on decay products
 745                     if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
 746                         && !KinematicSelection(iparticle, 1))
 747                     {
 748                         continue;
 749                     }
 750 //
 751 // Decay products
 752 // Select if mother was selected and is not tracked
 753
 754                     if (pSelected[ipa] &&
 755                         !trackIt[ipa]  &&     // mother will be  tracked ?
 756                         kfMo !=  5 &&         // mother is b-quark, don't store fragments
 757                         kfMo !=  4 &&         // mother is c-quark, don't store fragments
 758                         kf   != 92)           // don't store string
 759                     {
 760 //
 761 // Semi-stable or de-selected: diselect decay products:
 762 //
 763 //
 764                         if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
 765                         {
 766                             Int_t ipF = iparticle->GetFirstDaughter();
 767                             Int_t ipL = iparticle->GetLastDaughter();
 768                             if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
 769                         }
 770 //                      printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
 771                         pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
 772                     }
 773                 }
 774                 if (pSelected[i] == -1) pSelected[i] = 0;
 775                 if (!pSelected[i]) continue;
 776                 // Count quarks only if you did not include fragmentation
 777                 if (fFragmentation && kf <= 10) continue;
 778
 779                 nc++;
 780 // Decision on tracking
 781                 trackIt[i] = 0;
 782 //
 783 // Track final state particle
 784                 if (ks == 1) trackIt[i] = 1;
 785 // Track semi-stable particles
 786                 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
 787 // Track particles selected by process if undecayed.
 788                 if (fForceDecay == kNoDecay) {
 789                     if (ParentSelected(kf)) trackIt[i] = 1;
 790                 } else {
 791                     if (ParentSelected(kf)) trackIt[i] = 0;
 792                 }
 793                 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
 794 //
 795 //
 796
 797             } // particle selection loop
 798             if (nc > 0) {
 799                 for (i = 0; i<np; i++) {
 800                     if (!pSelected[i]) continue;
 801                     TParticle *  iparticle = (TParticle *) fParticles.At(i);
 802                     kf = CheckPDGCode(iparticle->GetPdgCode());
 803                     Int_t ks = iparticle->GetStatusCode();
 804                     p[0] = iparticle->Px();
 805                     p[1] = iparticle->Py();
 806                     p[2] = iparticle->Pz();
 807                     p[3] = iparticle->Energy();
 808
 809                     origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 810                     origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 811                     origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 812
 813                     Float_t tof   = kconv*iparticle->T();
 814                     Int_t ipa     = iparticle->GetFirstMother()-1;
 815                     Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
 816
 817                     PushTrack(fTrackIt*trackIt[i], iparent, kf,
 818                               p[0], p[1], p[2], p[3],
 819                               origin[0], origin[1], origin[2], tof,
 820                               polar[0], polar[1], polar[2],
 821                               kPPrimary, nt, 1., ks);
 822                     pParent[i] = nt;
 823                     KeepTrack(nt);
 824                     fNprimaries++;
 825                 } //  PushTrack loop
 826             }
 827         } else {
 828             nc = GenerateMB();
 829         } // mb ?
 830
 831         GetSubEventTime();
 832
 833         delete[] pParent;
 834         delete[] pSelected;
 835         delete[] trackIt;
 836
 837         if (nc > 0) {
 838           switch (fCountMode) {
 839           case kCountAll:
 840             // printf(" Count all \n");
 841             jev += nc;
 842             break;
 843           case kCountParents:
 844             // printf(" Count parents \n");
 845             jev += nParents;
 846             break;
 847           case kCountTrackables:
 848             // printf(" Count trackable \n");
 849             jev += nTkbles;
 850             break;
 851           }
 852             if (jev >= fNpart || fNpart == -1) {
 853                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 854
 855                 fQ  += fPythia->GetVINT(51);
 856                 fX1 += fPythia->GetVINT(41);
 857                 fX2 += fPythia->GetVINT(42);
 858                 fTrialsRun += fTrials;
 859                 fNev++;
 860                 MakeHeader();
 861                 break;
 862             }
 863         }
 864     } // event loop
 865     SetHighWaterMark(nt);
 866 //  adjust weight due to kinematic selection
 867 //    AdjustWeights();
 868 //  get cross-section
 869     fXsection=fPythia->GetPARI(1);
 870 }
 871
 872 Int_t  AliGenPythia::GenerateMB()
 873 {
 874 //
 875 // Min Bias selection and other global selections
 876 //
 877     Int_t i, kf, nt, iparent;
 878     Int_t nc = 0;
 879     Float_t p[4];
 880     Float_t polar[3]   =   {0,0,0};
 881     Float_t origin[3]  =   {0,0,0};
 882 //  converts from mm/c to s
 883     const Float_t kconv=0.001/2.999792458e8;
 884
 885
 886
 887     Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
 888
 889
 890
 891     Int_t* pParent = new Int_t[np];
 892     for (i=0; i< np; i++) pParent[i] = -1;
 893      if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
 894         TParticle* jet1 = (TParticle *) fParticles.At(6);
 895         TParticle* jet2 = (TParticle *) fParticles.At(7);
 896         if (!CheckTrigger(jet1, jet2)) {
 897           delete [] pParent;
 898           return 0;
 899         }
 900     }
 901
 902     // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
 903     if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
 904
 905       Bool_t ok = kFALSE;
 906
 907       Int_t pdg  = 0;
 908       if (fFragPhotonInCalo) pdg = 22   ; // Photon
 909       else if (fPi0InCalo) pdg = 111 ;    // Pi0
 910
 911       for (i=0; i< np; i++) {
 912         TParticle* iparticle = (TParticle *) fParticles.At(i);
 913         if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg &&
 914            iparticle->Pt() > fFragPhotonOrPi0MinPt){
 915           Int_t imother = iparticle->GetFirstMother() - 1;
 916           TParticle* pmother = (TParticle *) fParticles.At(imother);
 917           if(pdg == 111 ||
 918              (pdg == 22 && pmother->GetStatusCode() != 11)) //No photon from hadron decay
 919             {
 920               Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
 921               Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
 922               if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
 923                  (fCheckPHOS    && IsInPHOS(phi,eta)) )
 924                 ok =kTRUE;
 925             }
 926         }
 927       }
 928       if(!ok)
 929         return 0;
 930     }
 931
 932     // Select beauty jets with electron in EMCAL
 933     if (fProcess == kPyBeautyJets && fEleInEMCAL) {
 934
 935       Bool_t ok = kFALSE;
 936
 937       Int_t pdg  = 11; //electron
 938
 939       Float_t pt  = 0.;
 940       Float_t eta = 0.;
 941       Float_t phi = 0.;
 942       for (i=0; i< np; i++) {
 943         TParticle* iparticle = (TParticle *) fParticles.At(i);
 944         if(iparticle->GetStatusCode()==1 && TMath::Abs(iparticle->GetPdgCode())==pdg &&
 945            iparticle->Pt() > fElectronMinPt){
 946           pt = iparticle->Pt();
 947           phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
 948           eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
 949           if(IsInEMCAL(phi,eta))
 950             ok =kTRUE;
 951         }
 952       }
 953       if(!ok)
 954         return 0;
 955       AliDebug(5,Form("Found an electron jet (pt,eta,phi) = (%f,%f,%f)",pt,eta,phi));
 956     }
 957
 958     // Check for minimum multiplicity
 959     if (fTriggerMultiplicity > 0) {
 960       Int_t multiplicity = 0;
 961       for (i = 0; i < np; i++) {
 962         TParticle *  iparticle = (TParticle *) fParticles.At(i);
 963
 964         Int_t statusCode = iparticle->GetStatusCode();
 965
 966         // Initial state particle
 967         if (statusCode != 1)
 968           continue;
 969         // eta cut
 970         if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
 971           continue;
 972         // pt cut
 973         if (iparticle->Pt() < fTriggerMultiplicityPtMin)
 974             continue;
 975
 976         TParticlePDG* pdgPart = iparticle->GetPDG();
 977         if (pdgPart && pdgPart->Charge() == 0)
 978           continue;
 979
 980         ++multiplicity;
 981       }
 982
 983       if (multiplicity < fTriggerMultiplicity) {
 984         delete [] pParent;
 985         return 0;
 986       }
 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 == kPyMbWithDirectPhoton ||
1093           fProcess == kPyMbNonDiffr)
1094          && (fCutOnChild == 1) ) {
1095       if ( !CheckKinematicsOnChild() ) {
1096         delete[] pParent;
1097         return 0;
1098       }
1099     }
1100
1101
1102     for (i = 0; i < np; i++) {
1103         Int_t trackIt = 0;
1104         TParticle *  iparticle = (TParticle *) fParticles.At(i);
1105         kf = CheckPDGCode(iparticle->GetPdgCode());
1106         Int_t ks = iparticle->GetStatusCode();
1107         Int_t km = iparticle->GetFirstMother();
1108         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
1109             (ks != 1) ||
1110             ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)) {
1111             nc++;
1112             if (ks == 1) trackIt = 1;
1113             Int_t ipa = iparticle->GetFirstMother()-1;
1114
1115             iparent = (ipa > -1) ? pParent[ipa] : -1;
1116
1117 //
1118 // store track information
1119             p[0] = iparticle->Px();
1120             p[1] = iparticle->Py();
1121             p[2] = iparticle->Pz();
1122             p[3] = iparticle->Energy();
1123
1124
1125             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1126             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1127             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1128
1129             Float_t tof = fEventTime + kconv * iparticle->T();
1130
1131             PushTrack(fTrackIt*trackIt, iparent, kf,
1132                       p[0], p[1], p[2], p[3],
1133                       origin[0], origin[1], origin[2], tof,
1134                       polar[0], polar[1], polar[2],
1135                       kPPrimary, nt, 1., ks);
1136             fNprimaries++;
1137             KeepTrack(nt);
1138             pParent[i] = nt;
1139             SetHighWaterMark(nt);
1140
1141         } // select particle
1142     } // particle loop
1143
1144     delete[] pParent;
1145
1146     return 1;
1147 }
1148
1149
1150 void AliGenPythia::FinishRun()
1151 {
1152 // Print x-section summary
1153     fPythia->Pystat(1);
1154
1155     if (fNev > 0.) {
1156         fQ  /= fNev;
1157         fX1 /= fNev;
1158         fX2 /= fNev;
1159     }
1160
1161     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1162     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1163 }
1164
1165 void AliGenPythia::AdjustWeights() const
1166 {
1167 // Adjust the weights after generation of all events
1168 //
1169     if (gAlice) {
1170         TParticle *part;
1171         Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1172         for (Int_t i=0; i<ntrack; i++) {
1173             part= gAlice->GetMCApp()->Particle(i);
1174             part->SetWeight(part->GetWeight()*fKineBias);
1175         }
1176     }
1177 }
1178
1179 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
1180 {
1181 // Treat protons as inside nuclei with mass numbers a1 and a2
1182
1183     fAProjectile = a1;
1184     fATarget     = a2;
1185     fNucPdf      = pdfset;  // 0 EKS98 1 EPS08
1186     fSetNuclei   = kTRUE;
1187 }
1188
1189
1190 void AliGenPythia::MakeHeader()
1191 {
1192 //
1193 // Make header for the simulated event
1194 //
1195   if (gAlice) {
1196     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1197         gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1198   }
1199
1200 // Builds the event header, to be called after each event
1201     if (fHeader) delete fHeader;
1202     fHeader = new AliGenPythiaEventHeader("Pythia");
1203 //
1204 // Event type
1205     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1206 //
1207 // Number of trials
1208     ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1209 //
1210 // Event Vertex
1211     fHeader->SetPrimaryVertex(fVertex);
1212
1213 //
1214 // Number of primaries
1215     fHeader->SetNProduced(fNprimaries);
1216 //
1217 // Jets that have triggered
1218
1219     //Need to store jets for b-jet studies too!
1220     if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
1221     {
1222         Int_t ntrig, njet;
1223         Float_t jets[4][10];
1224         GetJets(njet, ntrig, jets);
1225
1226
1227         for (Int_t i = 0; i < ntrig; i++) {
1228             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1229                                                         jets[3][i]);
1230         }
1231     }
1232 //
1233 // Copy relevant information from external header, if present.
1234 //
1235     Float_t uqJet[4];
1236
1237     if (fRL) {
1238         AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1239         for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1240         {
1241             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
1242
1243
1244             exHeader->TriggerJet(i, uqJet);
1245             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1246         }
1247     }
1248 //
1249 // Store quenching parameters
1250 //
1251     if (fQuench){
1252         Double_t z[4];
1253         Double_t xp, yp;
1254         if (fQuench == 1) {
1255             // Pythia::Quench()
1256             fPythia->GetQuenchingParameters(xp, yp, z);
1257         } else if (fQuench == 2){
1258             // Pyquen
1259             Double_t r1 = PARIMP.rb1;
1260             Double_t r2 = PARIMP.rb2;
1261             Double_t b  = PARIMP.b1;
1262             Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1263             Double_t phi = PARIMP.psib1;
1264             xp = r * TMath::Cos(phi);
1265             yp = r * TMath::Sin(phi);
1266
1267         } else if (fQuench == 4) {
1268             // QPythia
1269             Double_t xy[2];
1270             Double_t i0i1[2];
1271             AliFastGlauber::Instance()->GetSavedXY(xy);
1272             AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
1273             xp = xy[0];
1274             yp = xy[1];
1275             ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
1276         }
1277
1278             ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1279             ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1280     }
1281 //
1282 // Store pt^hard
1283     ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1284 //
1285 //  Pass header
1286 //
1287     AddHeader(fHeader);
1288     fHeader = 0x0;
1289 }
1290
1291 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
1292 {
1293 // Check the kinematic trigger condition
1294 //
1295     Double_t eta[2];
1296     eta[0] = jet1->Eta();
1297     eta[1] = jet2->Eta();
1298     Double_t phi[2];
1299     phi[0] = jet1->Phi();
1300     phi[1] = jet2->Phi();
1301     Int_t    pdg[2];
1302     pdg[0] = jet1->GetPdgCode();
1303     pdg[1] = jet2->GetPdgCode();
1304     Bool_t   triggered = kFALSE;
1305
1306     if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess ==  kPyBeautyppMNRwmi) {
1307         Int_t njets = 0;
1308         Int_t ntrig = 0;
1309         Float_t jets[4][10];
1310 //
1311 // Use Pythia clustering on parton level to determine jet axis
1312 //
1313         GetJets(njets, ntrig, jets);
1314
1315         if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1316 //
1317     } else {
1318         Int_t ij = 0;
1319         Int_t ig = 1;
1320         if (pdg[0] == kGamma) {
1321             ij = 1;
1322             ig = 0;
1323         }
1324         //Check eta range first...
1325         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
1326             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1327         {
1328             //Eta is okay, now check phi range
1329             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
1330                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1331             {
1332                 triggered = kTRUE;
1333             }
1334         }
1335     }
1336     return triggered;
1337 }
1338
1339
1340
1341 Bool_t AliGenPythia::CheckKinematicsOnChild(){
1342 //
1343 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1344 //
1345     Bool_t checking = kFALSE;
1346     Int_t j, kcode, ks, km;
1347     Int_t nPartAcc = 0; //number of particles in the acceptance range
1348     Int_t numberOfAcceptedParticles = 1;
1349     if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1350     Int_t npart = fParticles.GetEntriesFast();
1351
1352     for (j = 0; j<npart; j++) {
1353         TParticle *  jparticle = (TParticle *) fParticles.At(j);
1354         kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1355         ks = jparticle->GetStatusCode();
1356         km = jparticle->GetFirstMother();
1357
1358         if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
1359             nPartAcc++;
1360         }
1361         if( numberOfAcceptedParticles <= nPartAcc){
1362           checking = kTRUE;
1363           break;
1364         }
1365     }
1366
1367     return checking;
1368 }
1369
1370 void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1371 {
1372   //
1373   // Load event into Pythia Common Block
1374   //
1375
1376   Int_t npart = stack -> GetNprimary();
1377   Int_t n0 = 0;
1378
1379   if (!flag) {
1380     (fPythia->GetPyjets())->N = npart;
1381   } else {
1382     n0 = (fPythia->GetPyjets())->N;
1383     (fPythia->GetPyjets())->N = n0 + npart;
1384   }
1385
1386
1387   for (Int_t part = 0; part < npart; part++) {
1388     TParticle *mPart = stack->Particle(part);
1389
1390     Int_t kf     =  mPart->GetPdgCode();
1391     Int_t ks     =  mPart->GetStatusCode();
1392     Int_t idf    =  mPart->GetFirstDaughter();
1393     Int_t idl    =  mPart->GetLastDaughter();
1394
1395     if (reHadr) {
1396             if (ks == 11 || ks == 12) {
1397         ks  -= 10;
1398         idf  = -1;
1399         idl  = -1;
1400             }
1401     }
1402
1403     Float_t px = mPart->Px();
1404     Float_t py = mPart->Py();
1405     Float_t pz = mPart->Pz();
1406     Float_t e  = mPart->Energy();
1407     Float_t m  = mPart->GetCalcMass();
1408
1409
1410     (fPythia->GetPyjets())->P[0][part+n0] = px;
1411     (fPythia->GetPyjets())->P[1][part+n0] = py;
1412     (fPythia->GetPyjets())->P[2][part+n0] = pz;
1413     (fPythia->GetPyjets())->P[3][part+n0] = e;
1414     (fPythia->GetPyjets())->P[4][part+n0] = m;
1415
1416     (fPythia->GetPyjets())->K[1][part+n0] = kf;
1417     (fPythia->GetPyjets())->K[0][part+n0] = ks;
1418     (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1419     (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1420     (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1421   }
1422 }
1423
1424 void  AliGenPythia::LoadEvent(TObjArray* stack, Int_t flag, Int_t reHadr)
1425 {
1426   //
1427   // Load event into Pythia Common Block
1428   //
1429
1430   Int_t npart = stack -> GetEntries();
1431   Int_t n0 = 0;
1432
1433   if (!flag) {
1434     (fPythia->GetPyjets())->N = npart;
1435   } else {
1436     n0 = (fPythia->GetPyjets())->N;
1437     (fPythia->GetPyjets())->N = n0 + npart;
1438   }
1439
1440
1441   for (Int_t part = 0; part < npart; part++) {
1442     TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
1443     Int_t kf     =  mPart->GetPdgCode();
1444     Int_t ks     =  mPart->GetStatusCode();
1445     Int_t idf    =  mPart->GetFirstDaughter();
1446     Int_t idl    =  mPart->GetLastDaughter();
1447
1448     if (reHadr) {
1449             if (ks == 11 || ks == 12) {
1450         ks  -= 10;
1451         idf  = -1;
1452         idl  = -1;
1453             }
1454     }
1455
1456     Float_t px = mPart->Px();
1457     Float_t py = mPart->Py();
1458     Float_t pz = mPart->Pz();
1459     Float_t e  = mPart->Energy();
1460     Float_t m  = mPart->GetCalcMass();
1461
1462
1463     (fPythia->GetPyjets())->P[0][part+n0] = px;
1464     (fPythia->GetPyjets())->P[1][part+n0] = py;
1465     (fPythia->GetPyjets())->P[2][part+n0] = pz;
1466     (fPythia->GetPyjets())->P[3][part+n0] = e;
1467     (fPythia->GetPyjets())->P[4][part+n0] = m;
1468
1469     (fPythia->GetPyjets())->K[1][part+n0] = kf;
1470     (fPythia->GetPyjets())->K[0][part+n0] = ks;
1471     (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1472     (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1473     (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1474   }
1475 }
1476
1477
1478 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1479 {
1480 //
1481 //  Calls the Pythia jet finding algorithm to find jets in the current event
1482 //
1483 //
1484 //
1485 //  Save jets
1486     Int_t n     = fPythia->GetN();
1487
1488 //
1489 //  Run Jet Finder
1490     fPythia->Pycell(njets);
1491     Int_t i;
1492     for (i = 0; i < njets; i++) {
1493         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1494         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1495         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1496         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1497
1498         jets[0][i] = px;
1499         jets[1][i] = py;
1500         jets[2][i] = pz;
1501         jets[3][i] = e;
1502     }
1503 }
1504
1505
1506
1507 void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1508 {
1509 //
1510 //  Calls the Pythia clustering algorithm to find jets in the current event
1511 //
1512     Int_t n     = fPythia->GetN();
1513     nJets       = 0;
1514     nJetsTrig   = 0;
1515     if (fJetReconstruction == kCluster) {
1516 //
1517 //  Configure cluster algorithm
1518 //
1519         fPythia->SetPARU(43, 2.);
1520         fPythia->SetMSTU(41, 1);
1521 //
1522 //  Call cluster algorithm
1523 //
1524         fPythia->Pyclus(nJets);
1525 //
1526 //  Loading jets from common block
1527 //
1528     } else {
1529
1530 //
1531 //  Run Jet Finder
1532         fPythia->Pycell(nJets);
1533     }
1534
1535     Int_t i;
1536     for (i = 0; i < nJets; i++) {
1537         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1538         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1539         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1540         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1541         Float_t pt    = TMath::Sqrt(px * px + py * py);
1542         Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);
1543         Float_t theta = TMath::ATan2(pt,pz);
1544         Float_t et    = e * TMath::Sin(theta);
1545         Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
1546         if (
1547             eta > fEtaMinJet && eta < fEtaMaxJet &&
1548             phi > fPhiMinJet && phi < fPhiMaxJet &&
1549             et  > fEtMinJet  && et  < fEtMaxJet
1550             )
1551         {
1552             jets[0][nJetsTrig] = px;
1553             jets[1][nJetsTrig] = py;
1554             jets[2][nJetsTrig] = pz;
1555             jets[3][nJetsTrig] = e;
1556             nJetsTrig++;
1557 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1558         } else {
1559 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1560         }
1561     }
1562 }
1563
1564 void AliGenPythia::GetSubEventTime()
1565 {
1566   // Calculates time of the next subevent
1567   fEventTime = 0.;
1568   if (fEventsTime) {
1569     TArrayF &array = *fEventsTime;
1570     fEventTime = array[fCurSubEvent++];
1571   }
1572   //  printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1573   return;
1574 }
1575
1576 Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta)
1577 {
1578   // Is particle in EMCAL acceptance?
1579   // phi in degrees, etamin=-etamax
1580   if(phi > fEMCALMinPhi  && phi < fEMCALMaxPhi &&
1581      eta < fEMCALEta  )
1582     return kTRUE;
1583   else
1584     return kFALSE;
1585 }
1586
1587 Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta)
1588 {
1589   // Is particle in PHOS acceptance?
1590   // Acceptance slightly larger considered.
1591   // phi in degrees, etamin=-etamax
1592   if(phi > fPHOSMinPhi  && phi < fPHOSMaxPhi &&
1593      eta < fPHOSEta  )
1594     return kTRUE;
1595   else
1596     return kFALSE;
1597 }
1598
1599 void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd)
1600 {
1601   //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1602   Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1603   Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1604   Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
1605
1606   //calculate deltaphi
1607   TParticle* ph = (TParticle *) fParticles.At(iphcand);
1608   Double_t phphi = ph->Phi();
1609   Double_t deltaphi = phiPHOS - phphi;
1610
1611
1612
1613   //loop for all particles and produce the phi rotation
1614   Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1615   Double_t oldphi, newphi;
1616   Double_t newVx, newVy, R, Vz, time;
1617   Double_t newPx, newPy, pt, Pz, e;
1618   for(Int_t i=0; i< np; i++) {
1619       TParticle* iparticle = (TParticle *) fParticles.At(i);
1620       oldphi = iparticle->Phi();
1621       newphi = oldphi + deltaphi;
1622       if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1623       if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1624
1625       R = iparticle->R();
1626       newVx = R*TMath::Cos(newphi);
1627       newVy = R*TMath::Sin(newphi);
1628       Vz = iparticle->Vz(); // don't transform
1629       time = iparticle->T(); // don't transform
1630
1631       pt = iparticle->Pt();
1632       newPx = pt*TMath::Cos(newphi);
1633       newPy = pt*TMath::Sin(newphi);
1634       Pz = iparticle->Pz(); // don't transform
1635       e = iparticle->Energy(); // don't transform
1636
1637       // apply rotation
1638       iparticle->SetProductionVertex(newVx, newVy, Vz, time);
1639       iparticle->SetMomentum(newPx, newPy, Pz, e);
1640
1641   } //end particle loop
1642
1643    // now let's check that we put correctly the candidate photon in PHOS
1644    Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1645    Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1646    if(IsInPHOS(phi,eta))
1647       okdd = kTRUE;
1648 }
1649
1650
1651 #ifdef never
1652 void AliGenPythia::Streamer(TBuffer &R__b)
1653 {
1654    // Stream an object of class AliGenPythia.
1655
1656    if (R__b.IsReading()) {
1657       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1658       AliGenerator::Streamer(R__b);
1659       R__b >> (Int_t&)fProcess;
1660       R__b >> (Int_t&)fStrucFunc;
1661       R__b >> (Int_t&)fForceDecay;
1662       R__b >> fEnergyCMS;
1663       R__b >> fKineBias;
1664       R__b >> fTrials;
1665       fParentSelect.Streamer(R__b);
1666       fChildSelect.Streamer(R__b);
1667       R__b >> fXsection;
1668 //      (AliPythia::Instance())->Streamer(R__b);
1669       R__b >> fPtHardMin;
1670       R__b >> fPtHardMax;
1671 //      if (fDecayer) fDecayer->Streamer(R__b);
1672    } else {
1673       R__b.WriteVersion(AliGenPythia::IsA());
1674       AliGenerator::Streamer(R__b);
1675       R__b << (Int_t)fProcess;
1676       R__b << (Int_t)fStrucFunc;
1677       R__b << (Int_t)fForceDecay;
1678       R__b << fEnergyCMS;
1679       R__b << fKineBias;
1680       R__b << fTrials;
1681       fParentSelect.Streamer(R__b);
1682       fChildSelect.Streamer(R__b);
1683       R__b << fXsection;
1684 //      R__b << fPythia;
1685       R__b << fPtHardMin;
1686       R__b << fPtHardMax;
1687       //     fDecayer->Streamer(R__b);
1688    }
1689 }
1690 #endif
1691
1692
1693