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     fCountMode(kCountAll),
 120     fHeader(0),
 121     fRL(0),
 122     fFileName(0),
 123     fFragPhotonInCalo(kFALSE),
 124     fPi0InCalo(kFALSE) ,
 125     fPhotonInCalo(kFALSE),
 126     fEleInEMCAL(kFALSE),
 127     fCheckEMCAL(kFALSE),
 128     fCheckPHOS(kFALSE),
 129     fCheckPHOSeta(kFALSE),
 130     fFragPhotonOrPi0MinPt(0),
 131     fPhotonMinPt(0),
 132     fElectronMinPt(0),
 133     fPHOSMinPhi(219.),
 134     fPHOSMaxPhi(321.),
 135     fPHOSEta(0.13),
 136     fEMCALMinPhi(79.),
 137     fEMCALMaxPhi(191.),
 138     fEMCALEta(0.71)
 139
 140 {
 141 // Default Constructor
 142   fEnergyCMS = 5500.;
 143   SetNuclei(0,0);
 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      fCountMode(kCountAll),
 216      fHeader(0),
 217      fRL(0),
 218      fFileName(0),
 219      fFragPhotonInCalo(kFALSE),
 220      fPi0InCalo(kFALSE) ,
 221      fPhotonInCalo(kFALSE),
 222      fEleInEMCAL(kFALSE),
 223      fCheckEMCAL(kFALSE),
 224      fCheckPHOS(kFALSE),
 225      fCheckPHOSeta(kFALSE),
 226      fFragPhotonOrPi0MinPt(0),
 227      fPhotonMinPt(0),
 228      fElectronMinPt(0),
 229      fPHOSMinPhi(219.),
 230      fPHOSMaxPhi(321.),
 231      fPHOSEta(0.13),
 232      fEMCALMinPhi(79.),
 233      fEMCALMaxPhi(191.),
 234      fEMCALEta(0.71)
 235 {
 236 // default charm production at 5. 5 TeV
 237 // semimuonic decay
 238 // structure function GRVHO
 239 //
 240     fEnergyCMS = 5500.;
 241     fName = "Pythia";
 242     fTitle= "Particle Generator using PYTHIA";
 243     SetForceDecay();
 244     // Set random number generator
 245     if (!AliPythiaRndm::GetPythiaRandom())
 246       AliPythiaRndm::SetPythiaRandom(GetRandom());
 247     SetNuclei(0,0);
 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         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     //
 894     //TO BE CHECKED:  Should we require this for Beauty Jets?
 895     //
 896     if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets) {
 897         TParticle* jet1 = (TParticle *) fParticles.At(6);
 898         TParticle* jet2 = (TParticle *) fParticles.At(7);
 899         if (!CheckTrigger(jet1, jet2)) {
 900           delete [] pParent;
 901           return 0;
 902         }
 903     }
 904
 905     // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
 906     if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
 907
 908       Bool_t ok = kFALSE;
 909
 910       Int_t pdg  = 0;
 911       if (fFragPhotonInCalo) pdg = 22   ; // Photon
 912       else if (fPi0InCalo) pdg = 111 ;    // Pi0
 913
 914       for (i=0; i< np; i++) {
 915         TParticle* iparticle = (TParticle *) fParticles.At(i);
 916         if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg &&
 917            iparticle->Pt() > fFragPhotonOrPi0MinPt){
 918           Int_t imother = iparticle->GetFirstMother() - 1;
 919           TParticle* pmother = (TParticle *) fParticles.At(imother);
 920           if(pdg == 111 ||
 921              (pdg == 22 && pmother->GetStatusCode() != 11)) //No photon from hadron decay
 922             {
 923               Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
 924               Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
 925               if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
 926                  (fCheckPHOS    && IsInPHOS(phi,eta)) )
 927                 ok =kTRUE;
 928             }
 929         }
 930       }
 931       if(!ok)
 932         return 0;
 933     }
 934
 935     // Select beauty jets with electron in EMCAL
 936     if (fProcess == kPyBeautyJets && fEleInEMCAL) {
 937
 938       Bool_t ok = kFALSE;
 939
 940       Int_t pdg  = 11; //electron
 941
 942       Float_t pt  = 0.;
 943       Float_t eta = 0.;
 944       Float_t phi = 0.;
 945       for (i=0; i< np; i++) {
 946         TParticle* iparticle = (TParticle *) fParticles.At(i);
 947         if(iparticle->GetStatusCode()==1 && TMath::Abs(iparticle->GetPdgCode())==pdg &&
 948            iparticle->Pt() > fElectronMinPt){
 949           pt = iparticle->Pt();
 950           phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
 951           eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
 952           if(IsInEMCAL(phi,eta))
 953             ok =kTRUE;
 954         }
 955       }
 956       if(!ok)
 957         return 0;
 958       AliDebug(5,Form("Found an electron jet (pt,eta,phi) = (%f,%f,%f)",pt,eta,phi));
 959     }
 960
 961     // Check for minimum multiplicity
 962     if (fTriggerMultiplicity > 0) {
 963       Int_t multiplicity = 0;
 964       for (i = 0; i < np; i++) {
 965         TParticle *  iparticle = (TParticle *) fParticles.At(i);
 966
 967         Int_t statusCode = iparticle->GetStatusCode();
 968
 969         // Initial state particle
 970         if (statusCode > 20)
 971           continue;
 972
 973         // skip quarks and gluons
 974         Int_t pdgCode = TMath::Abs(iparticle->GetPdgCode());
 975         if (pdgCode <= 10 || pdgCode == 21)
 976           continue;
 977
 978         if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
 979           continue;
 980
 981         TParticlePDG* pdgPart = iparticle->GetPDG();
 982         if (pdgPart && pdgPart->Charge() == 0)
 983           continue;
 984
 985         ++multiplicity;
 986       }
 987
 988       if (multiplicity < fTriggerMultiplicity) {
 989         delete [] pParent;
 990         return 0;
 991       }
 992
 993       Printf("Triggered on event with multiplicity of %d > %d", multiplicity, fTriggerMultiplicity);
 994     }
 995
 996      // Select events with a photon  pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
 997     if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
 998
 999       Bool_t okd = kFALSE;
1000
1001       Int_t pdg  = 22;
1002       Int_t iphcand = -1;
1003       for (i=0; i< np; i++) {
1004          TParticle* iparticle = (TParticle *) fParticles.At(i);
1005          Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
1006          Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
1007
1008          if(iparticle->GetStatusCode() == 1
1009             && iparticle->GetPdgCode() == pdg
1010             && iparticle->Pt() > fPhotonMinPt
1011             && eta < fPHOSEta){
1012
1013             // first check if the photon is in PHOS phi
1014             if(IsInPHOS(phi,eta)){
1015                 okd = kTRUE;
1016                 break;
1017             }
1018             if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi
1019
1020          }
1021       }
1022
1023       if(!okd && iphcand != -1) // execute rotation in phi
1024           RotatePhi(iphcand,okd);
1025
1026       if(!okd)
1027         return 0;
1028     }
1029
1030     if (fTriggerParticle) {
1031         Bool_t triggered = kFALSE;
1032         for (i = 0; i < np; i++) {
1033             TParticle *  iparticle = (TParticle *) fParticles.At(i);
1034             kf = CheckPDGCode(iparticle->GetPdgCode());
1035             if (kf != fTriggerParticle) continue;
1036             if (iparticle->Pt() == 0.) continue;
1037             if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
1038             triggered = kTRUE;
1039             break;
1040         }
1041         if (!triggered) {
1042           delete [] pParent;
1043           return 0;
1044         }
1045     }
1046
1047
1048     // Check if there is a ccbar or bbbar pair with at least one of the two
1049     // in fYMin < y < fYMax
1050     //
1051     // TO BE CHECKED:  Should we require this for beauty jets?
1052     //
1053     if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
1054       TParticle *partCheck;
1055       TParticle *mother;
1056       Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
1057       Bool_t  theChild=kFALSE;
1058       Float_t y;
1059       Int_t   pdg,mpdg,mpdgUpperFamily;
1060       for(i=0; i<np; i++) {
1061         partCheck = (TParticle*)fParticles.At(i);
1062         pdg = partCheck->GetPdgCode();
1063         if(TMath::Abs(pdg) == fFlavorSelect) { // quark
1064           if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1065           y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1066                              (partCheck->Energy()-partCheck->Pz()+1.e-13));
1067           if(y>fYMin && y<fYMax) inYcut=kTRUE;
1068         }
1069         if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1070           Int_t mi = partCheck->GetFirstMother() - 1;
1071           if(mi<0) continue;
1072           mother = (TParticle*)fParticles.At(mi);
1073           mpdg=TMath::Abs(mother->GetPdgCode());
1074           mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
1075           if ( ParentSelected(mpdg) ||
1076               (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1077             if (KinematicSelection(partCheck,1)) {
1078               theChild=kTRUE;
1079             }
1080           }
1081         }
1082       }
1083       if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
1084         delete[] pParent;
1085         return 0;
1086       }
1087       if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1088         delete[] pParent;
1089         return 0;
1090       }
1091
1092     }
1093
1094     //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
1095     if ( (fProcess == kPyW ||
1096           fProcess == kPyZ ||
1097           fProcess == kPyMbDefault ||
1098           fProcess == kPyMb ||
1099           fProcess == kPyMbWithDirectPhoton ||
1100           fProcess == kPyMbNonDiffr)
1101          && (fCutOnChild == 1) ) {
1102       if ( !CheckKinematicsOnChild() ) {
1103         delete[] pParent;
1104         return 0;
1105       }
1106     }
1107
1108
1109     for (i = 0; i < np; i++) {
1110         Int_t trackIt = 0;
1111         TParticle *  iparticle = (TParticle *) fParticles.At(i);
1112         kf = CheckPDGCode(iparticle->GetPdgCode());
1113         Int_t ks = iparticle->GetStatusCode();
1114         Int_t km = iparticle->GetFirstMother();
1115         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
1116             (ks != 1) ||
1117             //
1118             //TO BE CHECKED:  Should we require this for beauty jets?
1119             //
1120             ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)) {
1121             nc++;
1122             if (ks == 1) trackIt = 1;
1123             Int_t ipa = iparticle->GetFirstMother()-1;
1124
1125             iparent = (ipa > -1) ? pParent[ipa] : -1;
1126
1127 //
1128 // store track information
1129             p[0] = iparticle->Px();
1130             p[1] = iparticle->Py();
1131             p[2] = iparticle->Pz();
1132             p[3] = iparticle->Energy();
1133
1134
1135             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1136             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1137             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1138
1139             Float_t tof = fEventTime + kconv * iparticle->T();
1140
1141             PushTrack(fTrackIt*trackIt, iparent, kf,
1142                       p[0], p[1], p[2], p[3],
1143                       origin[0], origin[1], origin[2], tof,
1144                       polar[0], polar[1], polar[2],
1145                       kPPrimary, nt, 1., ks);
1146             fNprimaries++;
1147             KeepTrack(nt);
1148             pParent[i] = nt;
1149             SetHighWaterMark(nt);
1150
1151         } // select particle
1152     } // particle loop
1153
1154     delete[] pParent;
1155
1156     return 1;
1157 }
1158
1159
1160 void AliGenPythia::FinishRun()
1161 {
1162 // Print x-section summary
1163     fPythia->Pystat(1);
1164
1165     if (fNev > 0.) {
1166         fQ  /= fNev;
1167         fX1 /= fNev;
1168         fX2 /= fNev;
1169     }
1170
1171     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1172     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1173 }
1174
1175 void AliGenPythia::AdjustWeights() const
1176 {
1177 // Adjust the weights after generation of all events
1178 //
1179     if (gAlice) {
1180         TParticle *part;
1181         Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1182         for (Int_t i=0; i<ntrack; i++) {
1183             part= gAlice->GetMCApp()->Particle(i);
1184             part->SetWeight(part->GetWeight()*fKineBias);
1185         }
1186     }
1187 }
1188
1189 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
1190 {
1191 // Treat protons as inside nuclei with mass numbers a1 and a2
1192
1193     fAProjectile = a1;
1194     fATarget     = a2;
1195     fNucPdf      = pdfset;  // 0 EKS98 1 EPS08
1196     fSetNuclei   = kTRUE;
1197 }
1198
1199
1200 void AliGenPythia::MakeHeader()
1201 {
1202 //
1203 // Make header for the simulated event
1204 //
1205   if (gAlice) {
1206     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1207         gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1208   }
1209
1210 // Builds the event header, to be called after each event
1211     if (fHeader) delete fHeader;
1212     fHeader = new AliGenPythiaEventHeader("Pythia");
1213 //
1214 // Event type
1215     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1216 //
1217 // Number of trials
1218     ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1219 //
1220 // Event Vertex
1221     fHeader->SetPrimaryVertex(fVertex);
1222
1223 //
1224 // Number of primaries
1225     fHeader->SetNProduced(fNprimaries);
1226 //
1227 // Jets that have triggered
1228
1229     //Need to store jets for b-jet studies too!
1230     if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets)
1231     {
1232         Int_t ntrig, njet;
1233         Float_t jets[4][10];
1234         GetJets(njet, ntrig, jets);
1235
1236
1237         for (Int_t i = 0; i < ntrig; i++) {
1238             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1239                                                         jets[3][i]);
1240         }
1241     }
1242 //
1243 // Copy relevant information from external header, if present.
1244 //
1245     Float_t uqJet[4];
1246
1247     if (fRL) {
1248         AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1249         for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1250         {
1251             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
1252
1253
1254             exHeader->TriggerJet(i, uqJet);
1255             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1256         }
1257     }
1258 //
1259 // Store quenching parameters
1260 //
1261     if (fQuench){
1262         Double_t z[4];
1263         Double_t xp, yp;
1264         if (fQuench == 1) {
1265             // Pythia::Quench()
1266             fPythia->GetQuenchingParameters(xp, yp, z);
1267         } else if (fQuench == 2){
1268             // Pyquen
1269             Double_t r1 = PARIMP.rb1;
1270             Double_t r2 = PARIMP.rb2;
1271             Double_t b  = PARIMP.b1;
1272             Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1273             Double_t phi = PARIMP.psib1;
1274             xp = r * TMath::Cos(phi);
1275             yp = r * TMath::Sin(phi);
1276
1277         } else if (fQuench == 4) {
1278             // QPythia
1279             Double_t xy[2];
1280             Double_t i0i1[2];
1281             AliFastGlauber::Instance()->GetSavedXY(xy);
1282             AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
1283             xp = xy[0];
1284             yp = xy[1];
1285             ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
1286         }
1287
1288             ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1289             ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1290     }
1291 //
1292 // Store pt^hard
1293     ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1294 //
1295 //  Pass header
1296 //
1297     AddHeader(fHeader);
1298     fHeader = 0x0;
1299 }
1300
1301 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
1302 {
1303 // Check the kinematic trigger condition
1304 //
1305     Double_t eta[2];
1306     eta[0] = jet1->Eta();
1307     eta[1] = jet2->Eta();
1308     Double_t phi[2];
1309     phi[0] = jet1->Phi();
1310     phi[1] = jet2->Phi();
1311     Int_t    pdg[2];
1312     pdg[0] = jet1->GetPdgCode();
1313     pdg[1] = jet2->GetPdgCode();
1314     Bool_t   triggered = kFALSE;
1315
1316     //
1317     //TO BE CHECKED: If we call this method for kPyBeautyJets, we need it here
1318     //
1319     if (fProcess == kPyJets || fProcess == kPyBeautyJets) {
1320         Int_t njets = 0;
1321         Int_t ntrig = 0;
1322         Float_t jets[4][10];
1323 //
1324 // Use Pythia clustering on parton level to determine jet axis
1325 //
1326         GetJets(njets, ntrig, jets);
1327
1328         if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1329 //
1330     } else {
1331         Int_t ij = 0;
1332         Int_t ig = 1;
1333         if (pdg[0] == kGamma) {
1334             ij = 1;
1335             ig = 0;
1336         }
1337         //Check eta range first...
1338         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
1339             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1340         {
1341             //Eta is okay, now check phi range
1342             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
1343                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1344             {
1345                 triggered = kTRUE;
1346             }
1347         }
1348     }
1349     return triggered;
1350 }
1351
1352
1353
1354 Bool_t AliGenPythia::CheckKinematicsOnChild(){
1355 //
1356 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1357 //
1358     Bool_t checking = kFALSE;
1359     Int_t j, kcode, ks, km;
1360     Int_t nPartAcc = 0; //number of particles in the acceptance range
1361     Int_t numberOfAcceptedParticles = 1;
1362     if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1363     Int_t npart = fParticles.GetEntriesFast();
1364
1365     for (j = 0; j<npart; j++) {
1366         TParticle *  jparticle = (TParticle *) fParticles.At(j);
1367         kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1368         ks = jparticle->GetStatusCode();
1369         km = jparticle->GetFirstMother();
1370
1371         if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
1372             nPartAcc++;
1373         }
1374         if( numberOfAcceptedParticles <= nPartAcc){
1375           checking = kTRUE;
1376           break;
1377         }
1378     }
1379
1380     return checking;
1381 }
1382
1383 void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1384 {
1385   //
1386   // Load event into Pythia Common Block
1387   //
1388
1389   Int_t npart = stack -> GetNprimary();
1390   Int_t n0 = 0;
1391
1392   if (!flag) {
1393     (fPythia->GetPyjets())->N = npart;
1394   } else {
1395     n0 = (fPythia->GetPyjets())->N;
1396     (fPythia->GetPyjets())->N = n0 + npart;
1397   }
1398
1399
1400   for (Int_t part = 0; part < npart; part++) {
1401     TParticle *mPart = stack->Particle(part);
1402
1403     Int_t kf     =  mPart->GetPdgCode();
1404     Int_t ks     =  mPart->GetStatusCode();
1405     Int_t idf    =  mPart->GetFirstDaughter();
1406     Int_t idl    =  mPart->GetLastDaughter();
1407
1408     if (reHadr) {
1409             if (ks == 11 || ks == 12) {
1410         ks  -= 10;
1411         idf  = -1;
1412         idl  = -1;
1413             }
1414     }
1415
1416     Float_t px = mPart->Px();
1417     Float_t py = mPart->Py();
1418     Float_t pz = mPart->Pz();
1419     Float_t e  = mPart->Energy();
1420     Float_t m  = mPart->GetCalcMass();
1421
1422
1423     (fPythia->GetPyjets())->P[0][part+n0] = px;
1424     (fPythia->GetPyjets())->P[1][part+n0] = py;
1425     (fPythia->GetPyjets())->P[2][part+n0] = pz;
1426     (fPythia->GetPyjets())->P[3][part+n0] = e;
1427     (fPythia->GetPyjets())->P[4][part+n0] = m;
1428
1429     (fPythia->GetPyjets())->K[1][part+n0] = kf;
1430     (fPythia->GetPyjets())->K[0][part+n0] = ks;
1431     (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1432     (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1433     (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1434   }
1435 }
1436
1437 void  AliGenPythia::LoadEvent(TObjArray* stack, Int_t flag, Int_t reHadr)
1438 {
1439   //
1440   // Load event into Pythia Common Block
1441   //
1442
1443   Int_t npart = stack -> GetEntries();
1444   Int_t n0 = 0;
1445
1446   if (!flag) {
1447     (fPythia->GetPyjets())->N = npart;
1448   } else {
1449     n0 = (fPythia->GetPyjets())->N;
1450     (fPythia->GetPyjets())->N = n0 + npart;
1451   }
1452
1453
1454   for (Int_t part = 0; part < npart; part++) {
1455     TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
1456     Int_t kf     =  mPart->GetPdgCode();
1457     Int_t ks     =  mPart->GetStatusCode();
1458     Int_t idf    =  mPart->GetFirstDaughter();
1459     Int_t idl    =  mPart->GetLastDaughter();
1460
1461     if (reHadr) {
1462             if (ks == 11 || ks == 12) {
1463         ks  -= 10;
1464         idf  = -1;
1465         idl  = -1;
1466             }
1467     }
1468
1469     Float_t px = mPart->Px();
1470     Float_t py = mPart->Py();
1471     Float_t pz = mPart->Pz();
1472     Float_t e  = mPart->Energy();
1473     Float_t m  = mPart->GetCalcMass();
1474
1475
1476     (fPythia->GetPyjets())->P[0][part+n0] = px;
1477     (fPythia->GetPyjets())->P[1][part+n0] = py;
1478     (fPythia->GetPyjets())->P[2][part+n0] = pz;
1479     (fPythia->GetPyjets())->P[3][part+n0] = e;
1480     (fPythia->GetPyjets())->P[4][part+n0] = m;
1481
1482     (fPythia->GetPyjets())->K[1][part+n0] = kf;
1483     (fPythia->GetPyjets())->K[0][part+n0] = ks;
1484     (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1485     (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1486     (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1487   }
1488 }
1489
1490
1491 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1492 {
1493 //
1494 //  Calls the Pythia jet finding algorithm to find jets in the current event
1495 //
1496 //
1497 //
1498 //  Save jets
1499     Int_t n     = fPythia->GetN();
1500
1501 //
1502 //  Run Jet Finder
1503     fPythia->Pycell(njets);
1504     Int_t i;
1505     for (i = 0; i < njets; i++) {
1506         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1507         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1508         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1509         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1510
1511         jets[0][i] = px;
1512         jets[1][i] = py;
1513         jets[2][i] = pz;
1514         jets[3][i] = e;
1515     }
1516 }
1517
1518
1519
1520 void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1521 {
1522 //
1523 //  Calls the Pythia clustering algorithm to find jets in the current event
1524 //
1525     Int_t n     = fPythia->GetN();
1526     nJets       = 0;
1527     nJetsTrig   = 0;
1528     if (fJetReconstruction == kCluster) {
1529 //
1530 //  Configure cluster algorithm
1531 //
1532         fPythia->SetPARU(43, 2.);
1533         fPythia->SetMSTU(41, 1);
1534 //
1535 //  Call cluster algorithm
1536 //
1537         fPythia->Pyclus(nJets);
1538 //
1539 //  Loading jets from common block
1540 //
1541     } else {
1542
1543 //
1544 //  Run Jet Finder
1545         fPythia->Pycell(nJets);
1546     }
1547
1548     Int_t i;
1549     for (i = 0; i < nJets; i++) {
1550         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1551         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1552         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1553         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1554         Float_t pt    = TMath::Sqrt(px * px + py * py);
1555         Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);
1556         Float_t theta = TMath::ATan2(pt,pz);
1557         Float_t et    = e * TMath::Sin(theta);
1558         Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
1559         if (
1560             eta > fEtaMinJet && eta < fEtaMaxJet &&
1561             phi > fPhiMinJet && phi < fPhiMaxJet &&
1562             et  > fEtMinJet  && et  < fEtMaxJet
1563             )
1564         {
1565             jets[0][nJetsTrig] = px;
1566             jets[1][nJetsTrig] = py;
1567             jets[2][nJetsTrig] = pz;
1568             jets[3][nJetsTrig] = e;
1569             nJetsTrig++;
1570 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1571         } else {
1572 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1573         }
1574     }
1575 }
1576
1577 void AliGenPythia::GetSubEventTime()
1578 {
1579   // Calculates time of the next subevent
1580   fEventTime = 0.;
1581   if (fEventsTime) {
1582     TArrayF &array = *fEventsTime;
1583     fEventTime = array[fCurSubEvent++];
1584   }
1585   //  printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1586   return;
1587 }
1588
1589 Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta)
1590 {
1591   // Is particle in EMCAL acceptance?
1592   // phi in degrees, etamin=-etamax
1593   if(phi > fEMCALMinPhi  && phi < fEMCALMaxPhi &&
1594      eta < fEMCALEta  )
1595     return kTRUE;
1596   else
1597     return kFALSE;
1598 }
1599
1600 Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta)
1601 {
1602   // Is particle in PHOS acceptance?
1603   // Acceptance slightly larger considered.
1604   // phi in degrees, etamin=-etamax
1605   if(phi > fPHOSMinPhi  && phi < fPHOSMaxPhi &&
1606      eta < fPHOSEta  )
1607     return kTRUE;
1608   else
1609     return kFALSE;
1610 }
1611
1612 void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd)
1613 {
1614   //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1615   Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1616   Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1617   Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
1618
1619   //calculate deltaphi
1620   TParticle* ph = (TParticle *) fParticles.At(iphcand);
1621   Double_t phphi = ph->Phi();
1622   Double_t deltaphi = phiPHOS - phphi;
1623
1624
1625
1626   //loop for all particles and produce the phi rotation
1627   Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1628   Double_t oldphi, newphi;
1629   Double_t newVx, newVy, R, Vz, time;
1630   Double_t newPx, newPy, pt, Pz, e;
1631   for(Int_t i=0; i< np; i++) {
1632       TParticle* iparticle = (TParticle *) fParticles.At(i);
1633       oldphi = iparticle->Phi();
1634       newphi = oldphi + deltaphi;
1635       if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1636       if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1637
1638       R = iparticle->R();
1639       newVx = R*TMath::Cos(newphi);
1640       newVy = R*TMath::Sin(newphi);
1641       Vz = iparticle->Vz(); // don't transform
1642       time = iparticle->T(); // don't transform
1643
1644       pt = iparticle->Pt();
1645       newPx = pt*TMath::Cos(newphi);
1646       newPy = pt*TMath::Sin(newphi);
1647       Pz = iparticle->Pz(); // don't transform
1648       e = iparticle->Energy(); // don't transform
1649
1650       // apply rotation
1651       iparticle->SetProductionVertex(newVx, newVy, Vz, time);
1652       iparticle->SetMomentum(newPx, newPy, Pz, e);
1653
1654   } //end particle loop
1655
1656    // now let's check that we put correctly the candidate photon in PHOS
1657    Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1658    Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1659    if(IsInPHOS(phi,eta))
1660       okdd = kTRUE;
1661 }
1662
1663
1664 #ifdef never
1665 void AliGenPythia::Streamer(TBuffer &R__b)
1666 {
1667    // Stream an object of class AliGenPythia.
1668
1669    if (R__b.IsReading()) {
1670       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1671       AliGenerator::Streamer(R__b);
1672       R__b >> (Int_t&)fProcess;
1673       R__b >> (Int_t&)fStrucFunc;
1674       R__b >> (Int_t&)fForceDecay;
1675       R__b >> fEnergyCMS;
1676       R__b >> fKineBias;
1677       R__b >> fTrials;
1678       fParentSelect.Streamer(R__b);
1679       fChildSelect.Streamer(R__b);
1680       R__b >> fXsection;
1681 //      (AliPythia::Instance())->Streamer(R__b);
1682       R__b >> fPtHardMin;
1683       R__b >> fPtHardMax;
1684 //      if (fDecayer) fDecayer->Streamer(R__b);
1685    } else {
1686       R__b.WriteVersion(AliGenPythia::IsA());
1687       AliGenerator::Streamer(R__b);
1688       R__b << (Int_t)fProcess;
1689       R__b << (Int_t)fStrucFunc;
1690       R__b << (Int_t)fForceDecay;
1691       R__b << fEnergyCMS;
1692       R__b << fKineBias;
1693       R__b << fTrials;
1694       fParentSelect.Streamer(R__b);
1695       fChildSelect.Streamer(R__b);
1696       R__b << fXsection;
1697 //      R__b << fPythia;
1698       R__b << fPtHardMin;
1699       R__b << fPtHardMax;
1700       //     fDecayer->Streamer(R__b);
1701    }
1702 }
1703 #endif
1704
1705
1706