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