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         break;
 448     }
 449 //
 450 //
 451 //  JetFinder for Trigger
 452 //
 453 //  Configure detector (EMCAL like)
 454 //
 455     fPythia->SetPycellParameters(fPycellEtaMax,fPycellNEta, fPycellNPhi,
 456                                  fPycellThreshold, fPycellEtSeed,
 457                                  fPycellMinEtJet, fPycellMaxRadius);
 458 //
 459 //  This counts the total number of calls to Pyevnt() per run.
 460     fTrialsRun = 0;
 461     fQ         = 0.;
 462     fX1        = 0.;
 463     fX2        = 0.;
 464     fNev       = 0 ;
 465 //
 466 //
 467 //
 468     AliGenMC::Init();
 469 //
 470 //
 471 //
 472     if (fSetNuclei) {
 473         fDyBoost = 0;
 474         Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
 475     }
 476
 477     if (fQuench) {
 478         fPythia->InitQuenching(0., 0.1, 0.6e6, 0, 0.97, 30);
 479     }
 480
 481 //    fPythia->SetPARJ(200, 0.0);
 482
 483 //    if (fQuench == 3) {
 484 //      // Nestor's change of the splittings
 485 //      fPythia->SetPARJ(200, 0.8);
 486 //      fPythia->SetMSTJ(41, 1);  // QCD radiation only
 487 //      fPythia->SetMSTJ(42, 2);  // angular ordering
 488 //      fPythia->SetMSTJ(44, 2);  // option to run alpha_s
 489 //      fPythia->SetMSTJ(47, 0);  // No correction back to hard scattering element
 490 //      fPythia->SetMSTJ(50, 0);  // No coherence in first branching
 491 //      fPythia->SetPARJ(82, 1.); // Cut off for parton showers
 492 //    }
 493 }
 494
 495 void AliGenPythiaPlus::Generate()
 496 {
 497 // Generate one event
 498
 499     fDecayer->ForceDecay();
 500
 501     Double_t polar[3]   =   {0,0,0};
 502     Double_t origin[3]  =   {0,0,0};
 503     Double_t p[4];
 504 //  converts from mm/c to s
 505     const Float_t kconv=0.001/2.999792458e8;
 506 //
 507     Int_t nt=0;
 508     Int_t jev=0;
 509     Int_t j, kf;
 510     fTrials=0;
 511     fEventTime = 0.;
 512
 513
 514
 515     //  Set collision vertex position
 516     if (fVertexSmear == kPerEvent) Vertex();
 517
 518 //  event loop
 519     while(1)
 520     {
 521 //
 522 // Produce event
 523 //
 524 //
 525 // Switch hadronisation off
 526 //
 527 //      fPythia->SwitchHadronisationOff();
 528 //
 529 // Either produce new event or read partons from file
 530 //
 531         if (!fReadFromFile) {
 532             if (!fNewMIS) {
 533                 fPythia->GenerateEvent();
 534             } else {
 535                 fPythia->GenerateMIEvent();
 536             }
 537             fNpartons = fPythia->GetNumberOfParticles();
 538         } else {
 539             printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
 540             fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
 541             fPythia->SetNumberOfParticles(0);
 542             fPythia->LoadEvent(fRL->Stack(), 0 , 1);
 543             fPythia->EditEventList(21);
 544         }
 545
 546 //
 547 //  Run quenching routine
 548 //
 549         if (fQuench == 1) {
 550             fPythia->Quench();
 551         } else if (fQuench == 2){
 552             fPythia->Pyquen(208., 0, 0.);
 553         } else if (fQuench == 3) {
 554             // Quenching is via multiplicative correction of the splittings
 555         }
 556
 557 //
 558 // Switch hadronisation on
 559 //
 560 //      fPythia->SwitchHadronisationOn();
 561 //
 562 // .. and perform hadronisation
 563 //      printf("Calling hadronisation %d\n", fPythia->GetN());
 564 //      fPythia->HadronizeEvent();
 565         fTrials++;
 566         fPythia->GetParticles(&fParticles);
 567         Boost();
 568 //
 569 //
 570 //
 571         Int_t i;
 572
 573         fNprimaries = 0;
 574         Int_t np = fParticles.GetEntriesFast();
 575
 576         if (np == 0) continue;
 577 //
 578
 579 //
 580         Int_t* pParent   = new Int_t[np];
 581         Int_t* pSelected = new Int_t[np];
 582         Int_t* trackIt   = new Int_t[np];
 583         for (i = 0; i < np; i++) {
 584             pParent[i]   = -1;
 585             pSelected[i] =  0;
 586             trackIt[i]   =  0;
 587         }
 588
 589         Int_t nc = 0;        // Total n. of selected particles
 590         Int_t nParents = 0;  // Selected parents
 591         Int_t nTkbles = 0;   // Trackable particles
 592         if (fProcess != kPyMbDefault &&
 593             fProcess != kPyMb &&
 594             fProcess != kPyMbWithDirectPhoton &&
 595             fProcess != kPyJets &&
 596             fProcess != kPyDirectGamma &&
 597             fProcess != kPyMbNonDiffr  &&
 598             fProcess != kPyMbMSEL1     &&
 599             fProcess != kPyW &&
 600             fProcess != kPyZ &&
 601             fProcess != kPyCharmppMNRwmi &&
 602             fProcess != kPyBeautyppMNRwmi) {
 603
 604             for (i = 0; i < np; i++) {
 605                 TParticle* iparticle = (TParticle *) fParticles.At(i);
 606                 Int_t ks = iparticle->GetStatusCode();
 607                 kf = CheckPDGCode(iparticle->GetPdgCode());
 608 // No initial state partons
 609                 if (ks==21) continue;
 610 //
 611 // Heavy Flavor Selection
 612 //
 613                 // quark ?
 614                 kf = TMath::Abs(kf);
 615                 Int_t kfl = kf;
 616                 // Resonance
 617
 618                 if (kfl > 100000) kfl %= 100000;
 619                 if (kfl > 10000)  kfl %= 10000;
 620                 // meson ?
 621                 if  (kfl > 10) kfl/=100;
 622                 // baryon
 623                 if (kfl > 10) kfl/=10;
 624                 Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
 625                 Int_t kfMo = 0;
 626 //
 627 // Establish mother daughter relation between heavy quarks and mesons
 628 //
 629                 if (kf >= fFlavorSelect && kf <= 6) {
 630                     Int_t idau = (fPythia->Version() == 6) ? (iparticle->GetFirstDaughter() - 1) :(iparticle->GetFirstDaughter());
 631                     if (idau > -1) {
 632                         TParticle* daughter = (TParticle *) fParticles.At(idau);
 633                         Int_t pdgD = daughter->GetPdgCode();
 634                         if (pdgD == 91 || pdgD == 92) {
 635                             Int_t jmin = (fPythia->Version() == 6) ? (daughter->GetFirstDaughter() - 1) : (daughter->GetFirstDaughter());
 636                             Int_t jmax = (fPythia->Version() == 6) ? (daughter->GetLastDaughter() - 1)  : (daughter->GetLastDaughter());
 637
 638                             for (Int_t jp = jmin; jp <= jmax; jp++)
 639                                 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
 640                         } // is string or cluster
 641                     } // has daughter
 642                 } // heavy quark
 643
 644
 645                 if (ipa > -1) {
 646                     TParticle *  mother = (TParticle *) fParticles.At(ipa);
 647                     kfMo = TMath::Abs(mother->GetPdgCode());
 648                 }
 649
 650                 // What to keep in Stack?
 651                 Bool_t flavorOK = kFALSE;
 652                 Bool_t selectOK = kFALSE;
 653                 if (fFeedDownOpt) {
 654                     if (kfl >= fFlavorSelect) flavorOK = kTRUE;
 655                 } else {
 656                     if (kfl > fFlavorSelect) {
 657                         nc = -1;
 658                         break;
 659                     }
 660                     if (kfl == fFlavorSelect) flavorOK = kTRUE;
 661                 }
 662                 switch (fStackFillOpt) {
 663                 case kFlavorSelection:
 664                     selectOK = kTRUE;
 665                     break;
 666                 case kParentSelection:
 667                     if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
 668                     break;
 669                 }
 670                 if (flavorOK && selectOK) {
 671 //
 672 // Heavy flavor hadron or quark
 673 //
 674 // Kinematic seletion on final state heavy flavor mesons
 675                     if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
 676                     {
 677                         continue;
 678                     }
 679                     pSelected[i] = 1;
 680                     if (ParentSelected(kf)) ++nParents; // Update parent count
 681 //                  printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
 682                 } else {
 683 // Kinematic seletion on decay products
 684                     if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
 685                         && !KinematicSelection(iparticle, 1))
 686                     {
 687                         continue;
 688                     }
 689 //
 690 // Decay products
 691 // Select if mother was selected and is not tracked
 692
 693                     if (pSelected[ipa] &&
 694                         !trackIt[ipa]  &&     // mother will be  tracked ?
 695                         kfMo !=  5 &&         // mother is b-quark, don't store fragments
 696                         kfMo !=  4 &&         // mother is c-quark, don't store fragments
 697                         kf   != 92)           // don't store string
 698                     {
 699 //
 700 // Semi-stable or de-selected: diselect decay products:
 701 //
 702 //
 703                         if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
 704                         {
 705                             Int_t ipF = iparticle->GetFirstDaughter();
 706                             Int_t ipL = iparticle->GetLastDaughter();
 707                             if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
 708                         }
 709 //                      printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
 710                         pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
 711                     }
 712                 }
 713                 if (pSelected[i] == -1) pSelected[i] = 0;
 714                 if (!pSelected[i]) continue;
 715                 // Count quarks only if you did not include fragmentation
 716                 if (fFragmentation && kf <= 10) continue;
 717
 718                 nc++;
 719 // Decision on tracking
 720                 trackIt[i] = 0;
 721 //
 722 // Track final state particle
 723                 if (ks == 1) trackIt[i] = 1;
 724 // Track semi-stable particles
 725                 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
 726 // Track particles selected by process if undecayed.
 727                 if (fForceDecay == kNoDecay) {
 728                     if (ParentSelected(kf)) trackIt[i] = 1;
 729                 } else {
 730                     if (ParentSelected(kf)) trackIt[i] = 0;
 731                 }
 732                 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
 733 //
 734 //
 735
 736             } // particle selection loop
 737             if (nc > 0) {
 738                 for (i = 0; i < np; i++) {
 739                     if (!pSelected[i]) continue;
 740                     TParticle *  iparticle = (TParticle *) fParticles.At(i);
 741                     kf = CheckPDGCode(iparticle->GetPdgCode());
 742                     Int_t ks = iparticle->GetStatusCode();
 743                     p[0] = iparticle->Px();
 744                     p[1] = iparticle->Py();
 745                     p[2] = iparticle->Pz();
 746                     p[3] = iparticle->Energy();
 747
 748                     origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 749                     origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 750                     origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 751
 752                     Float_t tof   = fTime + kconv*iparticle->T();
 753                     Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
 754                     Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
 755
 756                     PushTrack(fTrackIt*trackIt[i], iparent, kf,
 757                               p[0], p[1], p[2], p[3],
 758                               origin[0], origin[1], origin[2], tof,
 759                               polar[0], polar[1], polar[2],
 760                               kPPrimary, nt, 1., ks);
 761                     pParent[i] = nt;
 762                     KeepTrack(nt);
 763                     fNprimaries++;
 764                 } //  PushTrack loop
 765             }
 766         } else {
 767             nc = GenerateMB();
 768         } // mb ?
 769
 770         GetSubEventTime();
 771
 772         delete[] pParent;
 773         delete[] pSelected;
 774         delete[] trackIt;
 775
 776         if (nc > 0) {
 777           switch (fCountMode) {
 778           case kCountAll:
 779             // printf(" Count all \n");
 780             jev += nc;
 781             break;
 782           case kCountParents:
 783             // printf(" Count parents \n");
 784             jev += nParents;
 785             break;
 786           case kCountTrackables:
 787             // printf(" Count trackable \n");
 788             jev += nTkbles;
 789             break;
 790           }
 791             if (jev >= fNpart || fNpart == -1) {
 792                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 793                 fPythia->GetXandQ(fQ, fX1, fX2);
 794                 fTrialsRun += fTrials;
 795                 fNev++;
 796                 MakeHeader();
 797                 break;
 798             }
 799         }
 800     } // event loop
 801     SetHighWaterMark(nt);
 802 //  Adjust weight due to kinematic selection
 803 //  AdjustWeights();
 804 //  get cross-section
 805     fXsection = fPythia->GetXSection();
 806 }
 807
 808 Int_t  AliGenPythiaPlus::GenerateMB()
 809 {
 810 //
 811 // Min Bias selection and other global selections
 812 //
 813     Int_t i, kf, nt, iparent;
 814     Int_t nc = 0;
 815     Float_t p[4];
 816     Float_t polar[3]   =   {0,0,0};
 817     Float_t origin[3]  =   {0,0,0};
 818 //  converts from mm/c to s
 819     const Float_t kconv = 0.001 / 2.999792458e8;
 820
 821     Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
 822
 823     Int_t* pParent = new Int_t[np];
 824     for (i=0; i< np; i++) pParent[i] = -1;
 825     if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
 826         TParticle* jet1 = (TParticle *) fParticles.At(6);
 827         TParticle* jet2 = (TParticle *) fParticles.At(7);
 828         if (!CheckTrigger(jet1, jet2)) {
 829           delete [] pParent;
 830           return 0;
 831         }
 832     }
 833
 834     // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
 835     if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
 836
 837       Bool_t ok = kFALSE;
 838
 839       Int_t pdg  = 0;
 840       if (fFragPhotonInCalo) pdg = 22   ; // Photon
 841       else if (fPi0InCalo) pdg = 111 ; // Pi0
 842
 843       for (i=0; i< np; i++) {
 844         TParticle* iparticle = (TParticle *) fParticles.At(i);
 845         if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg &&
 846            iparticle->Pt() > fFragPhotonOrPi0MinPt){
 847             Int_t imother = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
 848           TParticle* pmother = (TParticle *) fParticles.At(imother);
 849           if(pdg == 111 ||
 850              (pdg == 22 && pmother->GetStatusCode() != 11))//No photon from hadron decay
 851             {
 852               Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
 853               Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
 854               if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
 855                  (fCheckPHOS    && IsInPHOS(phi,eta)) )
 856                 ok =kTRUE;
 857             }
 858         }
 859       }
 860       if(!ok){
 861           delete [] pParent;
 862           return 0;
 863       }
 864     }
 865
 866
 867      // Select events with a photon  pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
 868     if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
 869
 870       Bool_t okd = kFALSE;
 871
 872       Int_t pdg  = 22;
 873       Int_t iphcand = -1;
 874       for (i=0; i< np; i++) {
 875          TParticle* iparticle = (TParticle *) fParticles.At(i);
 876          Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
 877          Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
 878
 879          if(iparticle->GetStatusCode() == 1
 880             && iparticle->GetPdgCode() == pdg
 881             && iparticle->Pt() > fPhotonMinPt
 882             && eta < fPHOSEta){
 883
 884             // first check if the photon is in PHOS phi
 885             if(IsInPHOS(phi,eta)){
 886                 okd = kTRUE;
 887                 break;
 888             }
 889             if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi
 890
 891          }
 892       }
 893
 894       if(!okd && iphcand != -1) // execute rotation in phi
 895           RotatePhi(iphcand,okd);
 896
 897       if(!okd) {
 898           delete[] pParent;
 899           return 0;
 900       }
 901     }
 902
 903     if (fTriggerParticle) {
 904         Bool_t triggered = kFALSE;
 905         for (i = 0; i < np; i++) {
 906             TParticle *  iparticle = (TParticle *) fParticles.At(i);
 907             kf = CheckPDGCode(iparticle->GetPdgCode());
 908             if (kf != fTriggerParticle) continue;
 909             if (iparticle->Pt() == 0.) continue;
 910             if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
 911             triggered = kTRUE;
 912             break;
 913         }
 914         if (!triggered) {
 915           delete [] pParent;
 916           return 0;
 917         }
 918     }
 919
 920
 921     // Check if there is a ccbar or bbbar pair with at least one of the two
 922     // in fYMin < y < fYMax
 923     if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
 924       TParticle *partCheck;
 925       TParticle *mother;
 926       Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
 927       Bool_t  theChild=kFALSE;
 928       Float_t y;
 929       Int_t   pdg,mpdg,mpdgUpperFamily;
 930       for(i=0; i<np; i++) {
 931         partCheck = (TParticle*)fParticles.At(i);
 932         pdg = partCheck->GetPdgCode();
 933         if(TMath::Abs(pdg) == fFlavorSelect) { // quark
 934           if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
 935           y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
 936                              (partCheck->Energy()-partCheck->Pz()+1.e-13));
 937           if(y>fYMin && y<fYMax) inYcut=kTRUE;
 938         }
 939
 940         if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
 941           Int_t mi = partCheck->GetFirstMother() - 1;
 942           if(mi<0) continue;
 943           mother = (TParticle*)fParticles.At(mi);
 944           mpdg=TMath::Abs(mother->GetPdgCode());
 945           mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
 946           if ( ParentSelected(mpdg) ||
 947               (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
 948             if (KinematicSelection(partCheck,1)) {
 949               theChild=kTRUE;
 950             }
 951           }
 952         }
 953       }
 954       if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
 955         delete[] pParent;
 956         return 0;
 957       }
 958       if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
 959         delete[] pParent;
 960         return 0;
 961       }
 962
 963     }
 964
 965     //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
 966     if ( (fProcess == kPyW ||
 967           fProcess == kPyZ ||
 968           fProcess == kPyMbDefault ||
 969           fProcess == kPyMb ||
 970           fProcess == kPyMbWithDirectPhoton ||
 971           fProcess == kPyMbNonDiffr)
 972          && (fCutOnChild == 1) ) {
 973       if ( !CheckKinematicsOnChild() ) {
 974         delete[] pParent;
 975         return 0;
 976       }
 977     }
 978
 979
 980     for (i = 0; i < np; i++) {
 981         Int_t trackIt = 0;
 982         TParticle *  iparticle = (TParticle *) fParticles.At(i);
 983         kf = CheckPDGCode(iparticle->GetPdgCode());
 984         Int_t ks = iparticle->GetStatusCode();
 985         Int_t km = iparticle->GetFirstMother();
 986         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
 987             (ks != 1) ||
 988             (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
 989             nc++;
 990             if (ks == 1) trackIt = 1;
 991
 992             Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
 993             iparent = (ipa > -1) ? pParent[ipa] : -1;
 994             if (ipa >= np) fPythia->EventListing();
 995
 996 //
 997 // store track information
 998             p[0] = iparticle->Px();
 999             p[1] = iparticle->Py();
1000             p[2] = iparticle->Pz();
1001             p[3] = iparticle->Energy();
1002
1003
1004             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1005             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1006             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1007
1008             Float_t tof = fTime + fEventTime + kconv * iparticle->T();
1009
1010             PushTrack(fTrackIt*trackIt, iparent, kf,
1011                       p[0], p[1], p[2], p[3],
1012                       origin[0], origin[1], origin[2], tof,
1013                       polar[0], polar[1], polar[2],
1014                       kPPrimary, nt, 1., ks);
1015             fNprimaries++;
1016
1017
1018             //
1019             // Special Treatment to store color-flow
1020             //
1021             /*
1022             if (ks == 3 || ks == 13 || ks == 14) {
1023                 TParticle* particle = 0;
1024                 if (fStack) {
1025                     particle = fStack->Particle(nt);
1026                 } else {
1027                     particle = AliRunLoader::Instance()->Stack()->Particle(nt);
1028                 }
1029                 particle->SetFirstDaughter(fPythia->GetK(2, i));
1030                 particle->SetLastDaughter(fPythia->GetK(3, i));
1031             }
1032             */
1033             KeepTrack(nt);
1034             pParent[i] = nt;
1035             SetHighWaterMark(nt);
1036
1037         } // select particle
1038     } // particle loop
1039
1040     delete[] pParent;
1041
1042     return 1;
1043 }
1044
1045
1046 void AliGenPythiaPlus::FinishRun()
1047 {
1048 // Print x-section summary
1049     fPythia->PrintStatistics();
1050
1051     if (fNev > 0.) {
1052         fQ  /= fNev;
1053         fX1 /= fNev;
1054         fX2 /= fNev;
1055     }
1056
1057     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1058     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1059 }
1060
1061 void AliGenPythiaPlus::AdjustWeights() const
1062 {
1063 // Adjust the weights after generation of all events
1064 //
1065     if (gAlice) {
1066         TParticle *part;
1067         Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1068         for (Int_t i=0; i<ntrack; i++) {
1069             part= gAlice->GetMCApp()->Particle(i);
1070             part->SetWeight(part->GetWeight()*fKineBias);
1071         }
1072     }
1073 }
1074
1075 void AliGenPythiaPlus::SetNuclei(Int_t a1, Int_t a2)
1076 {
1077 // Treat protons as inside nuclei with mass numbers a1 and a2
1078
1079     fAProjectile = a1;
1080     fATarget     = a2;
1081     fSetNuclei   = kTRUE;
1082 }
1083
1084
1085 void AliGenPythiaPlus::MakeHeader()
1086 {
1087 //
1088 // Make header for the simulated event
1089 //
1090   if (gAlice) {
1091     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1092         gAlice->GetEvNumber()<=fDebugEventLast) fPythia->EventListing();
1093   }
1094
1095 // Builds the event header, to be called after each event
1096     if (fHeader) delete fHeader;
1097     fHeader = new AliGenPythiaEventHeader("Pythia");
1098 //
1099 // Event type
1100     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->ProcessCode());
1101 //
1102 // Number of trials
1103     ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1104 //
1105 // Event Vertex
1106     fHeader->SetPrimaryVertex(fVertex);
1107     fHeader->SetInteractionTime(fTime+fEventTime);
1108 //
1109 // Number of primaries
1110     fHeader->SetNProduced(fNprimaries);
1111 //
1112 // Jets that have triggered
1113
1114     if (fProcess == kPyJets)
1115     {
1116         Int_t ntrig, njet;
1117         Float_t jets[4][10];
1118         GetJets(njet, ntrig, jets);
1119
1120
1121         for (Int_t i = 0; i < ntrig; i++) {
1122             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1123                                                         jets[3][i]);
1124         }
1125     }
1126 //
1127 // Copy relevant information from external header, if present.
1128 //
1129     Float_t uqJet[4];
1130
1131     if (fRL) {
1132         AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1133         for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1134         {
1135             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
1136
1137
1138             exHeader->TriggerJet(i, uqJet);
1139             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1140         }
1141     }
1142 //
1143 // Store quenching parameters
1144 //
1145     if (fQuench){
1146         Double_t z[4] = {0.};
1147         Double_t xp = 0.;
1148         Double_t yp = 0.;
1149         if (fQuench == 1) {
1150             // Pythia::Quench()
1151             fPythia->GetQuenchingParameters(xp, yp, z);
1152         } else {
1153             // Pyquen
1154             Double_t r1 = PARIMP.rb1;
1155             Double_t r2 = PARIMP.rb2;
1156             Double_t b  = PARIMP.b1;
1157             Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1158             Double_t phi = PARIMP.psib1;
1159             xp = r * TMath::Cos(phi);
1160             yp = r * TMath::Sin(phi);
1161
1162         }
1163             ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1164             ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1165         }
1166 //
1167 // Store pt^hard
1168     ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetPtHard());
1169 //
1170 //  Pass header
1171 //
1172     AddHeader(fHeader);
1173     fHeader = 0x0;
1174 }
1175
1176 Bool_t AliGenPythiaPlus::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
1177 {
1178 // Check the kinematic trigger condition
1179 //
1180     Double_t eta[2];
1181     eta[0] = jet1->Eta();
1182     eta[1] = jet2->Eta();
1183     Double_t phi[2];
1184     phi[0] = jet1->Phi();
1185     phi[1] = jet2->Phi();
1186     Int_t    pdg[2];
1187     pdg[0] = jet1->GetPdgCode();
1188     pdg[1] = jet2->GetPdgCode();
1189     Bool_t   triggered = kFALSE;
1190
1191     if (fProcess == kPyJets) {
1192         Int_t njets = 0;
1193         Int_t ntrig = 0;
1194         Float_t jets[4][10];
1195 //
1196 // Use Pythia clustering on parton level to determine jet axis
1197 //
1198         GetJets(njets, ntrig, jets);
1199
1200         if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1201 //
1202     } else {
1203         Int_t ij = 0;
1204         Int_t ig = 1;
1205         if (pdg[0] == kGamma) {
1206             ij = 1;
1207             ig = 0;
1208         }
1209         //Check eta range first...
1210         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
1211             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1212         {
1213             //Eta is okay, now check phi range
1214             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
1215                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1216             {
1217                 triggered = kTRUE;
1218             }
1219         }
1220     }
1221     return triggered;
1222 }
1223
1224
1225
1226 Bool_t AliGenPythiaPlus::CheckKinematicsOnChild(){
1227 //
1228 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1229 //
1230     Bool_t checking = kFALSE;
1231     Int_t j, kcode, ks, km;
1232     Int_t nPartAcc = 0; //number of particles in the acceptance range
1233     Int_t numberOfAcceptedParticles = 1;
1234     if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1235     Int_t npart = fParticles.GetEntriesFast();
1236
1237     for (j = 0; j<npart; j++) {
1238         TParticle *  jparticle = (TParticle *) fParticles.At(j);
1239         kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1240         ks = jparticle->GetStatusCode();
1241         km = jparticle->GetFirstMother();
1242
1243         if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
1244             nPartAcc++;
1245         }
1246         if( numberOfAcceptedParticles <= nPartAcc){
1247           checking = kTRUE;
1248           break;
1249         }
1250     }
1251
1252     return checking;
1253 }
1254
1255 void AliGenPythiaPlus::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1256 {
1257 //
1258 //  Calls the Pythia jet finding algorithm to find jets in the current event
1259 //
1260 //
1261 //
1262 //  Save jets
1263 //
1264 //  Run Jet Finder
1265     fPythia->Pycell(njets);
1266     Int_t i;
1267     for (i = 0; i < njets; i++) {
1268         Float_t px, py, pz, e;
1269         fPythia->GetJet(i, px, py, pz, e);
1270         jets[0][i] = px;
1271         jets[1][i] = py;
1272         jets[2][i] = pz;
1273         jets[3][i] = e;
1274     }
1275 }
1276
1277
1278 void  AliGenPythiaPlus::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1279 {
1280 //
1281 //  Calls the Pythia clustering algorithm to find jets in the current event
1282 //
1283     nJets       = 0;
1284     nJetsTrig   = 0;
1285     if (fJetReconstruction == kCluster) {
1286 //
1287 //  Configure cluster algorithm
1288 //
1289 //      fPythia->SetPARU(43, 2.);
1290 //      fPythia->SetMSTU(41, 1);
1291 //
1292 //  Call cluster algorithm
1293 //
1294         fPythia->Pyclus(nJets);
1295 //
1296 //  Loading jets from common block
1297 //
1298     } else {
1299
1300 //
1301 //  Run Jet Finder
1302         fPythia->Pycell(nJets);
1303     }
1304
1305     Int_t i;
1306     for (i = 0; i < nJets; i++) {
1307         Float_t px, py, pz, e;
1308         fPythia->GetJet(i, px, py, pz, e);
1309         Float_t pt    = TMath::Sqrt(px * px + py * py);
1310         Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);
1311         Float_t theta = TMath::ATan2(pt,pz);
1312         Float_t et    = e * TMath::Sin(theta);
1313         Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
1314         if (
1315             eta > fEtaMinJet && eta < fEtaMaxJet &&
1316             phi > fPhiMinJet && phi < fPhiMaxJet &&
1317             et  > fEtMinJet  && et  < fEtMaxJet
1318             )
1319         {
1320             jets[0][nJetsTrig] = px;
1321             jets[1][nJetsTrig] = py;
1322             jets[2][nJetsTrig] = pz;
1323             jets[3][nJetsTrig] = e;
1324             nJetsTrig++;
1325         } else {
1326         }
1327     }
1328 }
1329
1330 void AliGenPythiaPlus::GetSubEventTime()
1331 {
1332   // Calculates time of the next subevent
1333   fEventTime = 0.;
1334   if (fEventsTime) {
1335     TArrayF &array = *fEventsTime;
1336     fEventTime = array[fCurSubEvent++];
1337   }
1338   //  printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1339   return;
1340 }
1341
1342
1343
1344
1345 Bool_t AliGenPythiaPlus::IsInEMCAL(Float_t phi, Float_t eta) const
1346 {
1347   // Is particle in EMCAL acceptance?
1348   // phi in degrees, etamin=-etamax
1349   if(phi > fEMCALMinPhi  && phi < fEMCALMaxPhi &&
1350      eta < fEMCALEta  )
1351     return kTRUE;
1352   else
1353     return kFALSE;
1354 }
1355
1356 Bool_t AliGenPythiaPlus::IsInPHOS(Float_t phi, Float_t eta) const
1357 {
1358   // Is particle in PHOS acceptance?
1359   // Acceptance slightly larger considered.
1360   // phi in degrees, etamin=-etamax
1361   if(phi > fPHOSMinPhi  && phi < fPHOSMaxPhi &&
1362      eta < fPHOSEta  )
1363     return kTRUE;
1364   else
1365     return kFALSE;
1366 }
1367
1368 void AliGenPythiaPlus::RotatePhi(Int_t iphcand, Bool_t& okdd)
1369 {
1370   //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1371   Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1372   Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1373   Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
1374
1375   //calculate deltaphi
1376   TParticle* ph = (TParticle *) fParticles.At(iphcand);
1377   Double_t phphi = ph->Phi();
1378   Double_t deltaphi = phiPHOS - phphi;
1379
1380
1381
1382   //loop for all particles and produce the phi rotation
1383   Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1384   Double_t oldphi, newphi;
1385   Double_t newVx, newVy, R, Vz, time;
1386   Double_t newPx, newPy, pt, Pz, e;
1387   for(Int_t i=0; i< np; i++) {
1388       TParticle* iparticle = (TParticle *) fParticles.At(i);
1389       oldphi = iparticle->Phi();
1390       newphi = oldphi + deltaphi;
1391       if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1392       if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1393
1394       R = iparticle->R();
1395       newVx = R*TMath::Cos(newphi);
1396       newVy = R*TMath::Sin(newphi);
1397       Vz = iparticle->Vz(); // don't transform
1398       time = iparticle->T(); // don't transform
1399
1400       pt = iparticle->Pt();
1401       newPx = pt*TMath::Cos(newphi);
1402       newPy = pt*TMath::Sin(newphi);
1403       Pz = iparticle->Pz(); // don't transform
1404       e = iparticle->Energy(); // don't transform
1405
1406       // apply rotation
1407       iparticle->SetProductionVertex(newVx, newVy, Vz, time);
1408       iparticle->SetMomentum(newPx, newPy, Pz, e);
1409
1410   } //end particle loop
1411
1412    // now let's check that we put correctly the candidate photon in PHOS
1413    Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1414    Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1415    if(IsInPHOS(phi,eta))
1416       okdd = kTRUE;
1417 }
1418
1419
1420 #ifdef never
1421 void AliGenPythiaPlus::Streamer(TBuffer &R__b)
1422 {
1423    // Stream an object of class AliGenPythia.
1424
1425    if (R__b.IsReading()) {
1426       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1427       AliGenerator::Streamer(R__b);
1428       R__b >> (Int_t&)fProcess;
1429       R__b >> (Int_t&)fStrucFunc;
1430       R__b >> (Int_t&)fForceDecay;
1431       R__b >> fEnergyCMS;
1432       R__b >> fKineBias;
1433       R__b >> fTrials;
1434       fParentSelect.Streamer(R__b);
1435       fChildSelect.Streamer(R__b);
1436       R__b >> fXsection;
1437 //      (AliPythia::Instance())->Streamer(R__b);
1438       R__b >> fPtHardMin;
1439       R__b >> fPtHardMax;
1440 //      if (fDecayer) fDecayer->Streamer(R__b);
1441    } else {
1442       R__b.WriteVersion(AliGenPythiaPlus::IsA());
1443       AliGenerator::Streamer(R__b);
1444       R__b << (Int_t)fProcess;
1445       R__b << (Int_t)fStrucFunc;
1446       R__b << (Int_t)fForceDecay;
1447       R__b << fEnergyCMS;
1448       R__b << fKineBias;
1449       R__b << fTrials;
1450       fParentSelect.Streamer(R__b);
1451       fChildSelect.Streamer(R__b);
1452       R__b << fXsection;
1453 //      R__b << fPythia;
1454       R__b << fPtHardMin;
1455       R__b << fPtHardMax;
1456       //     fDecayer->Streamer(R__b);
1457    }
1458 }
1459 #endif
1460
1461
1462