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 <TDatabasePDG.h>
  29 #include <TParticle.h>
  30 #include <TPDGCode.h>
  31 #include <TSystem.h>
  32 #include <TTree.h>
  33 #include "AliConst.h"
  34 #include "AliDecayerPythia.h"
  35 #include "AliGenPythia.h"
  36 #include "AliHeader.h"
  37 #include "AliGenPythiaEventHeader.h"
  38 #include "AliPythia.h"
  39 #include "AliPythiaRndm.h"
  40 #include "AliRun.h"
  41 #include "AliStack.h"
  42 #include "AliRunLoader.h"
  43 #include "AliMC.h"
  44 #include "pyquenCommon.h"
  45
  46 ClassImp(AliGenPythia)
  47
  48 AliGenPythia::AliGenPythia()
  49                  :AliGenMC()
  50 {
  51 // Default Constructor
  52   fParticles = 0;
  53   fPythia    = 0;
  54   fHeader = 0;
  55   fReadFromFile = 0;
  56   fEventTime = 0.;
  57   fInteractionRate = 0.;
  58   fTimeWindow = 0.;
  59   fEventsTime = 0;
  60   fCurSubEvent = 0;
  61   fDecayer   = new AliDecayerPythia();
  62   SetEventListRange();
  63   SetJetPhiRange();
  64   SetJetEtaRange();
  65   SetJetEtRange();
  66   SetGammaPhiRange();
  67   SetGammaEtaRange();
  68   SetPtKick();
  69   SetQuench();
  70   SetHadronisation();
  71   SetTriggerParticle();
  72   SetNuclei(0,0);
  73   fSetNuclei = kFALSE;
  74   fNewMIS    = kFALSE;
  75   fHFoff     = kFALSE;
  76   fGinit     = 1;
  77   fGfinal    = 1;
  78
  79   if (!AliPythiaRndm::GetPythiaRandom())
  80     AliPythiaRndm::SetPythiaRandom(GetRandom());
  81 }
  82
  83 AliGenPythia::AliGenPythia(Int_t npart)
  84                  :AliGenMC(npart)
  85 {
  86 // default charm production at 5. 5 TeV
  87 // semimuonic decay
  88 // structure function GRVHO
  89 //
  90     fName = "Pythia";
  91     fTitle= "Particle Generator using PYTHIA";
  92     fXsection  = 0.;
  93     fReadFromFile = 0;
  94     fEventTime = 0.;
  95     fInteractionRate = 0.;
  96     fTimeWindow = 0.;
  97     fEventsTime = 0;
  98     fCurSubEvent = 0;
  99     SetProcess();
 100     SetStrucFunc();
 101     SetForceDecay();
 102     SetPtHard();
 103     SetYHard();
 104     SetEnergyCMS();
 105     fDecayer = new AliDecayerPythia();
 106     // Set random number generator
 107     if (!AliPythiaRndm::GetPythiaRandom())
 108       AliPythiaRndm::SetPythiaRandom(GetRandom());
 109     fFlavorSelect   = 0;
 110     // Produced particles
 111     fParticles = new TClonesArray("TParticle",1000);
 112     fHeader = 0;
 113     SetEventListRange();
 114     SetJetPhiRange();
 115     SetJetEtaRange();
 116     SetJetEtRange();
 117     SetGammaPhiRange();
 118     SetGammaEtaRange();
 119     SetJetReconstructionMode();
 120     SetQuench();
 121     SetHadronisation();
 122     SetPtKick();
 123     SetTriggerParticle();
 124     SetNuclei(0,0);
 125     // Options determining what to keep in the stack (Heavy flavour generation)
 126     fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
 127     fFeedDownOpt = kTRUE;             // allow feed down from higher family
 128     // Fragmentation on/off
 129     fFragmentation = kTRUE;
 130     // Default counting mode
 131     fCountMode = kCountAll;
 132     // Pycel
 133     SetPycellParameters();
 134     fSetNuclei = kFALSE;
 135     fNewMIS    = kFALSE;
 136     fHFoff     = kFALSE;
 137     fGinit     = 1;
 138     fGfinal    = 1;
 139  }
 140
 141 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
 142     :AliGenMC(Pythia)
 143 {
 144 // copy constructor
 145     Pythia.Copy(*this);
 146 }
 147
 148 AliGenPythia::~AliGenPythia()
 149 {
 150 // Destructor
 151   if(fEventsTime) delete fEventsTime;
 152 }
 153
 154 void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
 155 {
 156 // Generate pileup using user specified rate
 157     fInteractionRate = rate;
 158     fTimeWindow = timewindow;
 159     GeneratePileup();
 160 }
 161
 162 void AliGenPythia::GeneratePileup()
 163 {
 164 // Generate sub events time for pileup
 165     fEventsTime = 0;
 166     if(fInteractionRate == 0.) {
 167       Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
 168       return;
 169     }
 170
 171     Int_t npart = NumberParticles();
 172     if(npart < 0) {
 173       Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
 174       return;
 175     }
 176
 177     if(fEventsTime) delete fEventsTime;
 178     fEventsTime = new TArrayF(npart);
 179     TArrayF &array = *fEventsTime;
 180     for(Int_t ipart = 0; ipart < npart; ipart++)
 181       array[ipart] = 0.;
 182
 183     Float_t eventtime = 0.;
 184     while(1)
 185       {
 186         eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
 187         if(eventtime > fTimeWindow) break;
 188         array.Set(array.GetSize()+1);
 189         array[array.GetSize()-1] = eventtime;
 190       }
 191
 192     eventtime = 0.;
 193     while(1)
 194       {
 195         eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
 196         if(TMath::Abs(eventtime) > fTimeWindow) break;
 197         array.Set(array.GetSize()+1);
 198         array[array.GetSize()-1] = eventtime;
 199       }
 200
 201     SetNumberParticles(fEventsTime->GetSize());
 202 }
 203
 204 void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
 205                                        Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
 206 {
 207 // Set pycell parameters
 208     fPycellEtaMax    =  etamax;
 209     fPycellNEta      =  neta;
 210     fPycellNPhi      =  nphi;
 211     fPycellThreshold =  thresh;
 212     fPycellEtSeed    =  etseed;
 213     fPycellMinEtJet  =  minet;
 214     fPycellMaxRadius =  r;
 215 }
 216
 217
 218
 219 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
 220 {
 221   // Set a range of event numbers, for which a table
 222   // of generated particle will be printed
 223   fDebugEventFirst = eventFirst;
 224   fDebugEventLast  = eventLast;
 225   if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
 226 }
 227
 228 void AliGenPythia::Init()
 229 {
 230 // Initialisation
 231
 232     SetMC(AliPythia::Instance());
 233     fPythia=(AliPythia*) fMCEvGen;
 234
 235 //
 236     fParentWeight=1./Float_t(fNpart);
 237 //
 238 //  Forward Paramters to the AliPythia object
 239     fDecayer->SetForceDecay(fForceDecay);
 240     fDecayer->Init();
 241
 242
 243     fPythia->SetCKIN(3,fPtHardMin);
 244     fPythia->SetCKIN(4,fPtHardMax);
 245     fPythia->SetCKIN(7,fYHardMin);
 246     fPythia->SetCKIN(8,fYHardMax);
 247
 248     if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);
 249     // Fragmentation?
 250     if (fFragmentation) {
 251       fPythia->SetMSTP(111,1);
 252     } else {
 253       fPythia->SetMSTP(111,0);
 254     }
 255
 256
 257 //  initial state radiation
 258     fPythia->SetMSTP(61,fGinit);
 259 //  final state radiation
 260     fPythia->SetMSTP(71,fGfinal);
 261 //  pt - kick
 262     if (fPtKick > 0.) {
 263         fPythia->SetMSTP(91,1);
 264         fPythia->SetPARP(91,fPtKick);
 265     } else {
 266         fPythia->SetMSTP(91,0);
 267     }
 268
 269
 270     if (fReadFromFile) {
 271         fRL  =  AliRunLoader::Open(fFileName, "Partons");
 272         fRL->LoadKinematics();
 273         fRL->LoadHeader();
 274     } else {
 275         fRL = 0x0;
 276     }
 277 // Switch off Heavy Flavors on request
 278     if (fHFoff) {
 279         fPythia->SetMSTP(58, 3);
 280         fPythia->SetMSTJ(45, 3);
 281         for (Int_t i = 156; i <= 160; i++) fPythia->SetMDME(i, 1, 0);
 282     }
 283  //
 284     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
 285
 286 //  Parent and Children Selection
 287     switch (fProcess)
 288     {
 289     case kPyOldUEQ2ordered:
 290     case kPyOldUEQ2ordered2:
 291     case kPyOldPopcorn:
 292       break;
 293     case kPyCharm:
 294     case kPyCharmUnforced:
 295     case kPyCharmPbPbMNR:
 296     case kPyCharmpPbMNR:
 297     case kPyCharmppMNR:
 298     case kPyCharmppMNRwmi:
 299         fParentSelect[0] =   411;
 300         fParentSelect[1] =   421;
 301         fParentSelect[2] =   431;
 302         fParentSelect[3] =  4122;
 303         fFlavorSelect    =  4;
 304         break;
 305     case kPyD0PbPbMNR:
 306     case kPyD0pPbMNR:
 307     case kPyD0ppMNR:
 308         fParentSelect[0] =   421;
 309         fFlavorSelect    =   4;
 310         break;
 311     case kPyDPlusPbPbMNR:
 312     case kPyDPluspPbMNR:
 313     case kPyDPlusppMNR:
 314         fParentSelect[0] =   411;
 315         fFlavorSelect    =   4;
 316         break;
 317     case kPyDPlusStrangePbPbMNR:
 318     case kPyDPlusStrangepPbMNR:
 319     case kPyDPlusStrangeppMNR:
 320         fParentSelect[0] =   431;
 321         fFlavorSelect    =   4;
 322         break;
 323     case kPyBeauty:
 324     case kPyBeautyPbPbMNR:
 325     case kPyBeautypPbMNR:
 326     case kPyBeautyppMNR:
 327     case kPyBeautyppMNRwmi:
 328         fParentSelect[0]=  511;
 329         fParentSelect[1]=  521;
 330         fParentSelect[2]=  531;
 331         fParentSelect[3]= 5122;
 332         fParentSelect[4]= 5132;
 333         fParentSelect[5]= 5232;
 334         fParentSelect[6]= 5332;
 335         fFlavorSelect   = 5;
 336         break;
 337     case kPyBeautyUnforced:
 338         fParentSelect[0] =  511;
 339         fParentSelect[1] =  521;
 340         fParentSelect[2] =  531;
 341         fParentSelect[3] = 5122;
 342         fParentSelect[4] = 5132;
 343         fParentSelect[5] = 5232;
 344         fParentSelect[6] = 5332;
 345         fFlavorSelect    = 5;
 346         break;
 347     case kPyJpsiChi:
 348     case kPyJpsi:
 349         fParentSelect[0] = 443;
 350         break;
 351     case kPyMb:
 352     case kPyMbNonDiffr:
 353     case kPyMbMSEL1:
 354     case kPyJets:
 355     case kPyDirectGamma:
 356         break;
 357     case kPyW:
 358     case kPyZ:
 359         break;
 360     }
 361 //
 362 //
 363 //  JetFinder for Trigger
 364 //
 365 //  Configure detector (EMCAL like)
 366 //
 367     fPythia->SetPARU(51, fPycellEtaMax);
 368     fPythia->SetMSTU(51, fPycellNEta);
 369     fPythia->SetMSTU(52, fPycellNPhi);
 370 //
 371 //  Configure Jet Finder
 372 //
 373     fPythia->SetPARU(58,  fPycellThreshold);
 374     fPythia->SetPARU(52,  fPycellEtSeed);
 375     fPythia->SetPARU(53,  fPycellMinEtJet);
 376     fPythia->SetPARU(54,  fPycellMaxRadius);
 377     fPythia->SetMSTU(54,  2);
 378 //
 379 //  This counts the total number of calls to Pyevnt() per run.
 380     fTrialsRun = 0;
 381     fQ         = 0.;
 382     fX1        = 0.;
 383     fX2        = 0.;
 384     fNev       = 0 ;
 385 //
 386 //
 387 //
 388     AliGenMC::Init();
 389 //
 390 //
 391 //
 392     if (fSetNuclei) {
 393         fDyBoost = 0;
 394         Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
 395     }
 396
 397     if (fQuench) {
 398         fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
 399     }
 400 }
 401
 402 void AliGenPythia::Generate()
 403 {
 404 // Generate one event
 405
 406     fDecayer->ForceDecay();
 407
 408     Float_t polar[3]   =   {0,0,0};
 409     Float_t origin[3]  =   {0,0,0};
 410     Float_t p[4];
 411 //  converts from mm/c to s
 412     const Float_t kconv=0.001/2.999792458e8;
 413 //
 414     Int_t nt=0;
 415     Int_t jev=0;
 416     Int_t j, kf;
 417     fTrials=0;
 418     fEventTime = 0.;
 419
 420
 421
 422     //  Set collision vertex position
 423     if (fVertexSmear == kPerEvent) Vertex();
 424
 425 //  event loop
 426     while(1)
 427     {
 428 //
 429 // Produce event
 430 //
 431 //
 432 // Switch hadronisation off
 433 //
 434         fPythia->SetMSTJ(1, 0);
 435 //
 436 // Either produce new event or read partons from file
 437 //
 438         if (!fReadFromFile) {
 439             if (!fNewMIS) {
 440                 fPythia->Pyevnt();
 441             } else {
 442                 fPythia->Pyevnw();
 443             }
 444             fNpartons = fPythia->GetN();
 445         } else {
 446             printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
 447             fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
 448             fPythia->SetN(0);
 449             LoadEvent(fRL->Stack(), 0 , 1);
 450             fPythia->Pyedit(21);
 451         }
 452
 453 //
 454 //  Run quenching routine
 455 //
 456         if (fQuench == 1) {
 457             fPythia->Quench();
 458         } else if (fQuench == 2){
 459             fPythia->Pyquen(208., 0, 0.);
 460         }
 461 //
 462 // Switch hadronisation on
 463 //
 464         fPythia->SetMSTJ(1, 1);
 465 //
 466 // .. and perform hadronisation
 467 //      printf("Calling hadronisation %d\n", fPythia->GetN());
 468         fPythia->Pyexec();
 469         fTrials++;
 470         fPythia->ImportParticles(fParticles,"All");
 471         Boost();
 472 //
 473 //
 474 //
 475         Int_t i;
 476
 477
 478         Int_t np = fParticles->GetEntriesFast();
 479
 480         if (np == 0) continue;
 481 //
 482
 483 //
 484         Int_t* pParent   = new Int_t[np];
 485         Int_t* pSelected = new Int_t[np];
 486         Int_t* trackIt   = new Int_t[np];
 487         for (i = 0; i < np; i++) {
 488             pParent[i]   = -1;
 489             pSelected[i] =  0;
 490             trackIt[i]   =  0;
 491         }
 492
 493         Int_t nc = 0;        // Total n. of selected particles
 494         Int_t nParents = 0;  // Selected parents
 495         Int_t nTkbles = 0;   // Trackable particles
 496         if (fProcess != kPyMb && fProcess != kPyJets &&
 497             fProcess != kPyDirectGamma &&
 498             fProcess != kPyMbNonDiffr  &&
 499             fProcess != kPyMbMSEL1     &&
 500             fProcess != kPyW && fProcess != kPyZ &&
 501             fProcess != kPyCharmppMNRwmi && fProcess != kPyBeautyppMNRwmi) {
 502
 503             for (i = 0; i < np; i++) {
 504                 TParticle* iparticle = (TParticle *) fParticles->At(i);
 505                 Int_t ks = iparticle->GetStatusCode();
 506                 kf = CheckPDGCode(iparticle->GetPdgCode());
 507 // No initial state partons
 508                 if (ks==21) continue;
 509 //
 510 // Heavy Flavor Selection
 511 //
 512                 // quark ?
 513                 kf = TMath::Abs(kf);
 514                 Int_t kfl = kf;
 515                 // Resonance
 516
 517                 if (kfl > 100000) kfl %= 100000;
 518                 if (kfl > 10000)  kfl %= 10000;
 519                 // meson ?
 520                 if  (kfl > 10) kfl/=100;
 521                 // baryon
 522                 if (kfl > 10) kfl/=10;
 523                 Int_t ipa = iparticle->GetFirstMother()-1;
 524                 Int_t kfMo = 0;
 525 //
 526 // Establish mother daughter relation between heavy quarks and mesons
 527 //
 528                 if (kf >= fFlavorSelect && kf <= 6) {
 529                     Int_t idau = iparticle->GetFirstDaughter() - 1;
 530                     if (idau > -1) {
 531                         TParticle* daughter = (TParticle *) fParticles->At(idau);
 532                         Int_t pdgD = daughter->GetPdgCode();
 533                         if (pdgD == 91 || pdgD == 92) {
 534                             Int_t jmin = daughter->GetFirstDaughter() - 1;
 535                             Int_t jmax = daughter->GetLastDaughter()  - 1;
 536                             for (Int_t j = jmin; j <= jmax; j++)
 537                                 ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
 538                         } // is string or cluster
 539                     } // has daughter
 540                 } // heavy quark
 541
 542
 543                 if (ipa > -1) {
 544                     TParticle *  mother = (TParticle *) fParticles->At(ipa);
 545                     kfMo = TMath::Abs(mother->GetPdgCode());
 546                 }
 547
 548                 // What to keep in Stack?
 549                 Bool_t flavorOK = kFALSE;
 550                 Bool_t selectOK = kFALSE;
 551                 if (fFeedDownOpt) {
 552                     if (kfl >= fFlavorSelect) flavorOK = kTRUE;
 553                 } else {
 554                     if (kfl > fFlavorSelect) {
 555                         nc = -1;
 556                         break;
 557                     }
 558                     if (kfl == fFlavorSelect) flavorOK = kTRUE;
 559                 }
 560                 switch (fStackFillOpt) {
 561                 case kFlavorSelection:
 562                     selectOK = kTRUE;
 563                     break;
 564                 case kParentSelection:
 565                     if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
 566                     break;
 567                 }
 568                 if (flavorOK && selectOK) {
 569 //
 570 // Heavy flavor hadron or quark
 571 //
 572 // Kinematic seletion on final state heavy flavor mesons
 573                     if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
 574                     {
 575                         continue;
 576                     }
 577                     pSelected[i] = 1;
 578                     if (ParentSelected(kf)) ++nParents; // Update parent count
 579 //                  printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
 580                 } else {
 581 // Kinematic seletion on decay products
 582                     if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
 583                         && !KinematicSelection(iparticle, 1))
 584                     {
 585                         continue;
 586                     }
 587 //
 588 // Decay products
 589 // Select if mother was selected and is not tracked
 590
 591                     if (pSelected[ipa] &&
 592                         !trackIt[ipa]  &&     // mother will be  tracked ?
 593                         kfMo !=  5 &&         // mother is b-quark, don't store fragments
 594                         kfMo !=  4 &&         // mother is c-quark, don't store fragments
 595                         kf   != 92)           // don't store string
 596                     {
 597 //
 598 // Semi-stable or de-selected: diselect decay products:
 599 //
 600 //
 601                         if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
 602                         {
 603                             Int_t ipF = iparticle->GetFirstDaughter();
 604                             Int_t ipL = iparticle->GetLastDaughter();
 605                             if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
 606                         }
 607 //                      printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
 608                         pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
 609                     }
 610                 }
 611                 if (pSelected[i] == -1) pSelected[i] = 0;
 612                 if (!pSelected[i]) continue;
 613                 // Count quarks only if you did not include fragmentation
 614                 if (fFragmentation && kf <= 10) continue;
 615
 616                 nc++;
 617 // Decision on tracking
 618                 trackIt[i] = 0;
 619 //
 620 // Track final state particle
 621                 if (ks == 1) trackIt[i] = 1;
 622 // Track semi-stable particles
 623                 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
 624 // Track particles selected by process if undecayed.
 625                 if (fForceDecay == kNoDecay) {
 626                     if (ParentSelected(kf)) trackIt[i] = 1;
 627                 } else {
 628                     if (ParentSelected(kf)) trackIt[i] = 0;
 629                 }
 630                 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
 631 //
 632 //
 633
 634             } // particle selection loop
 635             if (nc > 0) {
 636                 for (i = 0; i<np; i++) {
 637                     if (!pSelected[i]) continue;
 638                     TParticle *  iparticle = (TParticle *) fParticles->At(i);
 639                     kf = CheckPDGCode(iparticle->GetPdgCode());
 640                     Int_t ks = iparticle->GetStatusCode();
 641                     p[0] = iparticle->Px();
 642                     p[1] = iparticle->Py();
 643                     p[2] = iparticle->Pz();
 644                     p[3] = iparticle->Energy();
 645
 646                     origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 647                     origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 648                     origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 649
 650                     Float_t tof   = kconv*iparticle->T();
 651                     Int_t ipa     = iparticle->GetFirstMother()-1;
 652                     Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
 653
 654                     PushTrack(fTrackIt*trackIt[i], iparent, kf,
 655                               p[0], p[1], p[2], p[3],
 656                               origin[0], origin[1], origin[2], tof,
 657                               polar[0], polar[1], polar[2],
 658                               kPPrimary, nt, 1., ks);
 659                     pParent[i] = nt;
 660                     KeepTrack(nt);
 661                 } //  PushTrack loop
 662             }
 663         } else {
 664             nc = GenerateMB();
 665         } // mb ?
 666
 667         GetSubEventTime();
 668
 669         delete[] pParent;
 670         delete[] pSelected;
 671         delete[] trackIt;
 672
 673         if (nc > 0) {
 674           switch (fCountMode) {
 675           case kCountAll:
 676             // printf(" Count all \n");
 677             jev += nc;
 678             break;
 679           case kCountParents:
 680             // printf(" Count parents \n");
 681             jev += nParents;
 682             break;
 683           case kCountTrackables:
 684             // printf(" Count trackable \n");
 685             jev += nTkbles;
 686             break;
 687           }
 688             if (jev >= fNpart || fNpart == -1) {
 689                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 690
 691                 fQ  += fPythia->GetVINT(51);
 692                 fX1 += fPythia->GetVINT(41);
 693                 fX2 += fPythia->GetVINT(42);
 694                 fTrialsRun += fTrials;
 695                 fNev++;
 696                 MakeHeader();
 697                 break;
 698             }
 699         }
 700     } // event loop
 701     SetHighWaterMark(nt);
 702 //  adjust weight due to kinematic selection
 703 //    AdjustWeights();
 704 //  get cross-section
 705     fXsection=fPythia->GetPARI(1);
 706 }
 707
 708 Int_t  AliGenPythia::GenerateMB()
 709 {
 710 //
 711 // Min Bias selection and other global selections
 712 //
 713     Int_t i, kf, nt, iparent;
 714     Int_t nc = 0;
 715     Float_t p[4];
 716     Float_t polar[3]   =   {0,0,0};
 717     Float_t origin[3]  =   {0,0,0};
 718 //  converts from mm/c to s
 719     const Float_t kconv=0.001/2.999792458e8;
 720
 721
 722
 723     Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
 724
 725
 726
 727     Int_t* pParent = new Int_t[np];
 728     for (i=0; i< np; i++) pParent[i] = -1;
 729     if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
 730         TParticle* jet1 = (TParticle *) fParticles->At(6);
 731         TParticle* jet2 = (TParticle *) fParticles->At(7);
 732         if (!CheckTrigger(jet1, jet2)) {
 733           delete [] pParent;
 734           return 0;
 735         }
 736     }
 737
 738     if (fTriggerParticle) {
 739         Bool_t triggered = kFALSE;
 740         for (i = 0; i < np; i++) {
 741             TParticle *  iparticle = (TParticle *) fParticles->At(i);
 742             kf = CheckPDGCode(iparticle->GetPdgCode());
 743             if (kf != fTriggerParticle) continue;
 744             if (iparticle->Pt() == 0.) continue;
 745             if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
 746             triggered = kTRUE;
 747             break;
 748         }
 749         if (!triggered) {
 750           delete [] pParent;
 751           return 0;
 752         }
 753     }
 754
 755
 756     // Check if there is a ccbar or bbbar pair with at least one of the two
 757     // in fYMin < y < fYMax
 758     if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
 759       TParticle *hvq;
 760       Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
 761       Float_t yQ;
 762       Int_t   pdgQ;
 763       for(i=0; i<np; i++) {
 764         hvq = (TParticle*)fParticles->At(i);
 765         pdgQ = hvq->GetPdgCode();
 766         if(TMath::Abs(pdgQ) != fFlavorSelect) continue;
 767         if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
 768         yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/
 769                             (hvq->Energy()-hvq->Pz()+1.e-13));
 770         if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE;
 771       }
 772       if (!theQ || !theQbar || !inYcut) {
 773         delete[] pParent;
 774         return 0;
 775       }
 776     }
 777
 778     //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
 779     if ( (fProcess == kPyW || fProcess == kPyZ || fProcess == kPyMb || fProcess == kPyMbNonDiffr)
 780          && (fCutOnChild == 1) ) {
 781       if ( !CheckKinematicsOnChild() ) {
 782         delete[] pParent;
 783         return 0;
 784       }
 785     }
 786
 787
 788     for (i = 0; i < np; i++) {
 789         Int_t trackIt = 0;
 790         TParticle *  iparticle = (TParticle *) fParticles->At(i);
 791         kf = CheckPDGCode(iparticle->GetPdgCode());
 792         Int_t ks = iparticle->GetStatusCode();
 793         Int_t km = iparticle->GetFirstMother();
 794         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
 795             (ks != 1) ||
 796             (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
 797             nc++;
 798             if (ks == 1) trackIt = 1;
 799             Int_t ipa = iparticle->GetFirstMother()-1;
 800
 801             iparent = (ipa > -1) ? pParent[ipa] : -1;
 802
 803 //
 804 // store track information
 805             p[0] = iparticle->Px();
 806             p[1] = iparticle->Py();
 807             p[2] = iparticle->Pz();
 808             p[3] = iparticle->Energy();
 809
 810
 811             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 812             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 813             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 814
 815             Float_t tof = fEventTime + kconv * iparticle->T();
 816
 817             PushTrack(fTrackIt*trackIt, iparent, kf,
 818                       p[0], p[1], p[2], p[3],
 819                       origin[0], origin[1], origin[2], tof,
 820                       polar[0], polar[1], polar[2],
 821                       kPPrimary, nt, 1., ks);
 822             //
 823             // Special Treatment to store color-flow
 824             //
 825             if (ks == 3 || ks == 13 || ks == 14) {
 826                 TParticle* particle = 0;
 827                 if (fStack) {
 828                     particle = fStack->Particle(nt);
 829                 } else {
 830                     particle = gAlice->Stack()->Particle(nt);
 831                 }
 832                 particle->SetFirstDaughter(fPythia->GetK(2, i));
 833                 particle->SetLastDaughter(fPythia->GetK(3, i));
 834             }
 835
 836             KeepTrack(nt);
 837             pParent[i] = nt;
 838             SetHighWaterMark(nt);
 839
 840         } // select particle
 841     } // particle loop
 842
 843     delete[] pParent;
 844
 845     return 1;
 846 }
 847
 848
 849 void AliGenPythia::FinishRun()
 850 {
 851 // Print x-section summary
 852     fPythia->Pystat(1);
 853
 854     if (fNev > 0.) {
 855         fQ  /= fNev;
 856         fX1 /= fNev;
 857         fX2 /= fNev;
 858     }
 859
 860     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
 861     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
 862 }
 863
 864 void AliGenPythia::AdjustWeights() const
 865 {
 866 // Adjust the weights after generation of all events
 867 //
 868     if (gAlice) {
 869         TParticle *part;
 870         Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
 871         for (Int_t i=0; i<ntrack; i++) {
 872             part= gAlice->GetMCApp()->Particle(i);
 873             part->SetWeight(part->GetWeight()*fKineBias);
 874         }
 875     }
 876 }
 877
 878 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
 879 {
 880 // Treat protons as inside nuclei with mass numbers a1 and a2
 881
 882     fAProjectile = a1;
 883     fATarget     = a2;
 884     fSetNuclei   = kTRUE;
 885 }
 886
 887
 888 void AliGenPythia::MakeHeader()
 889 {
 890 //
 891 // Make header for the simulated event
 892 //
 893   if (gAlice) {
 894     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
 895         gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
 896   }
 897
 898 // Builds the event header, to be called after each event
 899     if (fHeader) delete fHeader;
 900     fHeader = new AliGenPythiaEventHeader("Pythia");
 901 //
 902 // Event type
 903     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
 904 //
 905 // Number of trials
 906     ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
 907 //
 908 // Event Vertex
 909     fHeader->SetPrimaryVertex(fVertex);
 910 //
 911 // Jets that have triggered
 912
 913     if (fProcess == kPyJets)
 914     {
 915         Int_t ntrig, njet;
 916         Float_t jets[4][10];
 917         GetJets(njet, ntrig, jets);
 918
 919
 920         for (Int_t i = 0; i < ntrig; i++) {
 921             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
 922                                                         jets[3][i]);
 923         }
 924     }
 925 //
 926 // Copy relevant information from external header, if present.
 927 //
 928     Float_t uqJet[4];
 929
 930     if (fRL) {
 931         AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
 932         for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
 933         {
 934             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
 935
 936
 937             exHeader->TriggerJet(i, uqJet);
 938             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
 939         }
 940     }
 941 //
 942 // Store quenching parameters
 943 //
 944     if (fQuench){
 945         Double_t z[4];
 946         Double_t xp, yp;
 947         if (fQuench == 1) {
 948             // Pythia::Quench()
 949             fPythia->GetQuenchingParameters(xp, yp, z);
 950         } else {
 951             // Pyquen
 952             Double_t r1 = PARIMP.rb1;
 953             Double_t r2 = PARIMP.rb2;
 954             Double_t b  = PARIMP.b1;
 955             Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
 956             Double_t phi = PARIMP.psib1;
 957             xp = r * TMath::Cos(phi);
 958             yp = r * TMath::Sin(phi);
 959
 960         }
 961             ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
 962             ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
 963         }
 964 //
 965 // Store pt^hard
 966     ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
 967 //
 968 //  Pass header
 969 //
 970     AddHeader(fHeader);
 971 }
 972
 973 void AliGenPythia::AddHeader(AliGenEventHeader* header)
 974 {
 975     // Add header to container or runloader
 976     if (fContainer) {
 977         fContainer->AddHeader(header);
 978     } else {
 979         AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
 980     }
 981 }
 982
 983
 984 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
 985 {
 986 // Check the kinematic trigger condition
 987 //
 988     Double_t eta[2];
 989     eta[0] = jet1->Eta();
 990     eta[1] = jet2->Eta();
 991     Double_t phi[2];
 992     phi[0] = jet1->Phi();
 993     phi[1] = jet2->Phi();
 994     Int_t    pdg[2];
 995     pdg[0] = jet1->GetPdgCode();
 996     pdg[1] = jet2->GetPdgCode();
 997     Bool_t   triggered = kFALSE;
 998
 999     if (fProcess == kPyJets) {
1000         Int_t njets = 0;
1001         Int_t ntrig = 0;
1002         Float_t jets[4][10];
1003 //
1004 // Use Pythia clustering on parton level to determine jet axis
1005 //
1006         GetJets(njets, ntrig, jets);
1007
1008         if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1009 //
1010     } else {
1011         Int_t ij = 0;
1012         Int_t ig = 1;
1013         if (pdg[0] == kGamma) {
1014             ij = 1;
1015             ig = 0;
1016         }
1017         //Check eta range first...
1018         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
1019             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1020         {
1021             //Eta is okay, now check phi range
1022             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
1023                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1024             {
1025                 triggered = kTRUE;
1026             }
1027         }
1028     }
1029     return triggered;
1030 }
1031
1032
1033
1034 Bool_t AliGenPythia::CheckKinematicsOnChild(){
1035 //
1036 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1037 //
1038     Bool_t checking = kFALSE;
1039     Int_t j, kcode, ks, km;
1040     Int_t nPartAcc = 0; //number of particles in the acceptance range
1041     Int_t numberOfAcceptedParticles = 1;
1042     if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1043     Int_t npart = fParticles->GetEntriesFast();
1044
1045     for (j = 0; j<npart; j++) {
1046         TParticle *  jparticle = (TParticle *) fParticles->At(j);
1047         kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1048         ks = jparticle->GetStatusCode();
1049         km = jparticle->GetFirstMother();
1050
1051         if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
1052             nPartAcc++;
1053         }
1054         if( numberOfAcceptedParticles <= nPartAcc){
1055           checking = kTRUE;
1056           break;
1057         }
1058     }
1059
1060     return checking;
1061 }
1062
1063
1064 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
1065 {
1066 // Assignment operator
1067     rhs.Copy(*this);
1068     return *this;
1069 }
1070
1071 void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1072 {
1073 //
1074 // Load event into Pythia Common Block
1075 //
1076
1077     Int_t npart = stack -> GetNprimary();
1078     Int_t n0 = 0;
1079
1080     if (!flag) {
1081         (fPythia->GetPyjets())->N = npart;
1082     } else {
1083         n0 = (fPythia->GetPyjets())->N;
1084         (fPythia->GetPyjets())->N = n0 + npart;
1085     }
1086
1087
1088     for (Int_t part = 0; part < npart; part++) {
1089         TParticle *mPart = stack->Particle(part);
1090
1091         Int_t kf     =  mPart->GetPdgCode();
1092         Int_t ks     =  mPart->GetStatusCode();
1093         Int_t idf    =  mPart->GetFirstDaughter();
1094         Int_t idl    =  mPart->GetLastDaughter();
1095
1096         if (reHadr) {
1097             if (ks == 11 || ks == 12) {
1098                 ks  -= 10;
1099                 idf  = -1;
1100                 idl  = -1;
1101             }
1102         }
1103
1104         Float_t px = mPart->Px();
1105         Float_t py = mPart->Py();
1106         Float_t pz = mPart->Pz();
1107         Float_t e  = mPart->Energy();
1108         Float_t m  = mPart->GetCalcMass();
1109
1110
1111         (fPythia->GetPyjets())->P[0][part+n0] = px;
1112         (fPythia->GetPyjets())->P[1][part+n0] = py;
1113         (fPythia->GetPyjets())->P[2][part+n0] = pz;
1114         (fPythia->GetPyjets())->P[3][part+n0] = e;
1115         (fPythia->GetPyjets())->P[4][part+n0] = m;
1116
1117         (fPythia->GetPyjets())->K[1][part+n0] = kf;
1118         (fPythia->GetPyjets())->K[0][part+n0] = ks;
1119         (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1120         (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1121         (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1122     }
1123 }
1124
1125
1126 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1127 {
1128 //
1129 //  Calls the Pythia jet finding algorithm to find jets in the current event
1130 //
1131 //
1132 //
1133 //  Save jets
1134     Int_t n     = fPythia->GetN();
1135
1136 //
1137 //  Run Jet Finder
1138     fPythia->Pycell(njets);
1139     Int_t i;
1140     for (i = 0; i < njets; i++) {
1141         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1142         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1143         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1144         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1145
1146         jets[0][i] = px;
1147         jets[1][i] = py;
1148         jets[2][i] = pz;
1149         jets[3][i] = e;
1150     }
1151 }
1152
1153
1154
1155 void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1156 {
1157 //
1158 //  Calls the Pythia clustering algorithm to find jets in the current event
1159 //
1160     Int_t n     = fPythia->GetN();
1161     nJets       = 0;
1162     nJetsTrig   = 0;
1163     if (fJetReconstruction == kCluster) {
1164 //
1165 //  Configure cluster algorithm
1166 //
1167         fPythia->SetPARU(43, 2.);
1168         fPythia->SetMSTU(41, 1);
1169 //
1170 //  Call cluster algorithm
1171 //
1172         fPythia->Pyclus(nJets);
1173 //
1174 //  Loading jets from common block
1175 //
1176     } else {
1177
1178 //
1179 //  Run Jet Finder
1180         fPythia->Pycell(nJets);
1181     }
1182
1183     Int_t i;
1184     for (i = 0; i < nJets; i++) {
1185         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1186         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1187         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1188         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1189         Float_t pt    = TMath::Sqrt(px * px + py * py);
1190         Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);
1191         Float_t theta = TMath::ATan2(pt,pz);
1192         Float_t et    = e * TMath::Sin(theta);
1193         Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
1194         if (
1195             eta > fEtaMinJet && eta < fEtaMaxJet &&
1196             phi > fPhiMinJet && phi < fPhiMaxJet &&
1197             et  > fEtMinJet  && et  < fEtMaxJet
1198             )
1199         {
1200             jets[0][nJetsTrig] = px;
1201             jets[1][nJetsTrig] = py;
1202             jets[2][nJetsTrig] = pz;
1203             jets[3][nJetsTrig] = e;
1204             nJetsTrig++;
1205 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1206         } else {
1207 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1208         }
1209     }
1210 }
1211
1212 void AliGenPythia::GetSubEventTime()
1213 {
1214   // Calculates time of the next subevent
1215   fEventTime = 0.;
1216   if (fEventsTime) {
1217     TArrayF &array = *fEventsTime;
1218     fEventTime = array[fCurSubEvent++];
1219   }
1220   //  printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1221   return;
1222 }
1223
1224 #ifdef never
1225 void AliGenPythia::Streamer(TBuffer &R__b)
1226 {
1227    // Stream an object of class AliGenPythia.
1228
1229    if (R__b.IsReading()) {
1230       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1231       AliGenerator::Streamer(R__b);
1232       R__b >> (Int_t&)fProcess;
1233       R__b >> (Int_t&)fStrucFunc;
1234       R__b >> (Int_t&)fForceDecay;
1235       R__b >> fEnergyCMS;
1236       R__b >> fKineBias;
1237       R__b >> fTrials;
1238       fParentSelect.Streamer(R__b);
1239       fChildSelect.Streamer(R__b);
1240       R__b >> fXsection;
1241 //      (AliPythia::Instance())->Streamer(R__b);
1242       R__b >> fPtHardMin;
1243       R__b >> fPtHardMax;
1244 //      if (fDecayer) fDecayer->Streamer(R__b);
1245    } else {
1246       R__b.WriteVersion(AliGenPythia::IsA());
1247       AliGenerator::Streamer(R__b);
1248       R__b << (Int_t)fProcess;
1249       R__b << (Int_t)fStrucFunc;
1250       R__b << (Int_t)fForceDecay;
1251       R__b << fEnergyCMS;
1252       R__b << fKineBias;
1253       R__b << fTrials;
1254       fParentSelect.Streamer(R__b);
1255       fChildSelect.Streamer(R__b);
1256       R__b << fXsection;
1257 //      R__b << fPythia;
1258       R__b << fPtHardMin;
1259       R__b << fPtHardMax;
1260       //     fDecayer->Streamer(R__b);
1261    }
1262 }
1263 #endif
1264
1265
1266