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