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