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