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