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