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