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