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