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