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