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