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 /*
  17 $Log$
  18 Revision 1.2  2003/03/27 09:28:38  morsch
  19 Charm production process configuration for pPb and pp added
  20 (A. Dainese, N. Carrer)
  21
  22 Revision 1.1  2003/03/15 15:00:48  morsch
  23 Classed imported from EVGEN.
  24
  25 Revision 1.69  2003/01/14 10:50:19  alibrary
  26 Cleanup of STEER coding conventions
  27
  28 Revision 1.68  2002/12/11 09:16:16  morsch
  29 Use GetJets to fill header.
  30
  31 Revision 1.67  2002/12/09 15:24:09  morsch
  32 Same trigger routine can use Pycell or Pyclus.
  33
  34 Revision 1.66  2002/12/09 08:22:56  morsch
  35 UA1 jet finder (Pycell) for software triggering added.
  36
  37 Revision 1.65  2002/11/19 08:57:10  morsch
  38 Configuration of pt-kick added.
  39
  40 Revision 1.64  2002/11/15 00:43:06  morsch
  41 Changes for kPyJets
  42 - initial and final state g-radiation + pt-kick default
  43 - trigger based on parton clusters (using pyclus)
  44 - trigger jets are stored in header.
  45
  46 Revision 1.63  2002/10/14 14:55:35  hristov
  47 Merging the VirtualMC branch to the main development branch (HEAD)
  48
  49 Revision 1.52.4.4  2002/10/10 16:40:08  hristov
  50 Updating VirtualMC to v3-09-02
  51
  52 Revision 1.62  2002/09/24 10:00:01  morsch
  53 CheckTrigger() corrected.
  54
  55 Revision 1.61  2002/07/30 09:52:38  morsch
  56 Call SetGammaPhiRange() and SetGammaEtaRange() in the constructor.
  57
  58 Revision 1.60  2002/07/19 14:49:03  morsch
  59 Typo corrected.
  60
  61 Revision 1.59  2002/07/19 14:35:36  morsch
  62 Count total number of trials. Print mean Q, x1, x2.
  63
  64 Revision 1.58  2002/07/17 10:04:09  morsch
  65 SetYHard method added.
  66
  67 Revision 1.57  2002/05/22 13:22:53  morsch
  68 Process kPyMbNonDiffr added.
  69
  70 Revision 1.56  2002/04/26 10:30:01  morsch
  71 Option kPyBeautyPbMNR added. (N. Carrer).
  72
  73 Revision 1.55  2002/04/17 10:23:56  morsch
  74 Coding Rule violations corrected.
  75
  76 Revision 1.54  2002/03/28 11:49:10  morsch
  77 Pass status code in SetTrack.
  78
  79 Revision 1.53  2002/03/25 14:51:13  morsch
  80 New stack-fill and count options introduced (N. Carrer).
  81
  82 Revision 1.51  2002/03/06 08:46:57  morsch
  83 - Loop until np-1
  84 - delete dyn. alloc. arrays (N. Carrer)
  85
  86 Revision 1.50  2002/03/03 13:48:50  morsch
  87 Option  kPyCharmPbMNR added. Produce charm pairs in agreement with MNR
  88 NLO calculations (Nicola Carrer).
  89
  90 Revision 1.49  2002/02/08 16:50:50  morsch
  91 Add name and title in constructor.
  92
  93 Revision 1.48  2001/12/20 11:44:28  morsch
  94 Add kinematic bias for direct gamma production.
  95
  96 Revision 1.47  2001/12/19 14:45:00  morsch
  97 Store number of trials in header.
  98
  99 Revision 1.46  2001/12/19 10:36:19  morsch
 100 Add possibility if jet kinematic biasing.
 101
 102 Revision 1.45  2001/11/28 08:06:52  morsch
 103 Use fMaxLifeTime parameter.
 104
 105 Revision 1.44  2001/11/27 13:13:07  morsch
 106 Maximum lifetime for long-lived particles to be put on the stack is parameter.
 107 It can be set via SetMaximumLifetime(..).
 108
 109 Revision 1.43  2001/10/21 18:35:56  hristov
 110 Several pointers were set to zero in the default constructors to avoid memory management problems
 111
 112 Revision 1.42  2001/10/15 08:21:55  morsch
 113 Vertex truncation settings moved to AliGenMC.
 114
 115 Revision 1.41  2001/10/08 08:45:42  morsch
 116 Possibility of vertex cut added.
 117
 118 Revision 1.40  2001/09/25 11:30:23  morsch
 119 Pass event vertex to header.
 120
 121 Revision 1.39  2001/07/27 17:09:36  morsch
 122 Use local SetTrack, KeepTrack and SetHighWaterMark methods
 123 to delegate either to local stack or to stack owned by AliRun.
 124 (Piotr Skowronski, A.M.)
 125
 126 Revision 1.38  2001/07/13 10:58:54  morsch
 127 - Some coded moved to AliGenMC
 128 - Improved handling of secondary vertices.
 129
 130 Revision 1.37  2001/06/28 11:17:28  morsch
 131 SetEventListRange setter added. Events in specified range are listed for
 132 debugging. (Yuri Kharlov)
 133
 134 Revision 1.36  2001/03/30 07:05:49  morsch
 135 Final print-out in finish run.
 136 Write parton system for jet-production (preliminary solution).
 137
 138 Revision 1.35  2001/03/09 13:03:40  morsch
 139 Process_t and Struc_Func_t moved to AliPythia.h
 140
 141 Revision 1.34  2001/02/14 15:50:40  hristov
 142 The last particle in event marked using SetHighWaterMark
 143
 144 Revision 1.33  2001/01/30 09:23:12  hristov
 145 Streamers removed (R.Brun)
 146
 147 Revision 1.32  2001/01/26 19:55:51  hristov
 148 Major upgrade of AliRoot code
 149
 150 Revision 1.31  2001/01/17 10:54:31  hristov
 151 Better protection against FPE
 152
 153 Revision 1.30  2000/12/18 08:55:35  morsch
 154 Make AliPythia dependent generartors work with new scheme of random number generation
 155
 156 Revision 1.29  2000/12/04 11:22:03  morsch
 157 Init of sRandom as in 1.15
 158
 159 Revision 1.28  2000/12/02 11:41:39  morsch
 160 Use SetRandom() to initialize random number generator in constructor.
 161
 162 Revision 1.27  2000/11/30 20:29:02  morsch
 163 Initialise static variable sRandom in constructor: sRandom = fRandom;
 164
 165 Revision 1.26  2000/11/30 07:12:50  alibrary
 166 Introducing new Rndm and QA classes
 167
 168 Revision 1.25  2000/10/18 19:11:27  hristov
 169 Division by zero fixed
 170
 171 Revision 1.24  2000/09/18 10:41:35  morsch
 172 Add possibility to use nuclear structure functions from PDF library V8.
 173
 174 Revision 1.23  2000/09/14 14:05:40  morsch
 175 dito
 176
 177 Revision 1.22  2000/09/14 14:02:22  morsch
 178 - Correct conversion from mm to cm when passing particle vertex to MC.
 179 - Correct handling of fForceDecay == all.
 180
 181 Revision 1.21  2000/09/12 14:14:55  morsch
 182 Call fDecayer->ForceDecay() at the beginning of Generate().
 183
 184 Revision 1.20  2000/09/06 14:29:33  morsch
 185 Use AliPythia for event generation an AliDecayPythia for decays.
 186 Correct handling of "nodecay" option
 187
 188 Revision 1.19  2000/07/11 18:24:56  fca
 189 Coding convention corrections + few minor bug fixes
 190
 191 Revision 1.18  2000/06/30 12:40:34  morsch
 192 Pythia takes care of vertex smearing. Correct conversion from Pythia units (mm) to
 193 Geant units (cm).
 194
 195 Revision 1.17  2000/06/09 20:34:07  morsch
 196 All coding rule violations except RS3 corrected
 197
 198 Revision 1.16  2000/05/15 15:04:20  morsch
 199 The full event is written for fNtrack = -1
 200 Coding rule violations corrected.
 201
 202 Revision 1.15  2000/04/26 10:14:24  morsch
 203 Particles array has one entry more than pythia particle list. Upper bound of
 204 particle loop changed to np-1 (R. Guernane, AM)
 205
 206 Revision 1.14  2000/04/05 08:36:13  morsch
 207 Check status code of particles in Pythia event
 208 to avoid double counting as partonic state and final state particle.
 209
 210 Revision 1.13  1999/11/09 07:38:48  fca
 211 Changes for compatibility with version 2.23 of ROOT
 212
 213 Revision 1.12  1999/11/03 17:43:20  fca
 214 New version from G.Martinez & A.Morsch
 215
 216 Revision 1.11  1999/09/29 09:24:14  fca
 217 Introduction of the Copyright and cvs Log
 218 */
 219
 220 //
 221 // Generator using the TPythia interface (via AliPythia)
 222 // to generate pp collisions.
 223 // Using SetNuclei() also nuclear modifications to the structure functions
 224 // can be taken into account. This makes, of course, only sense for the
 225 // generation of the products of hard processes (heavy flavor, jets ...)
 226 //
 227 // andreas.morsch@cern.ch
 228 //
 229
 230 #include <TDatabasePDG.h>
 231 #include <TParticle.h>
 232 #include <TPDGCode.h>
 233 #include <TSystem.h>
 234 #include <TTree.h>
 235
 236 #include "AliConst.h"
 237 #include "AliDecayerPythia.h"
 238 #include "AliGenPythia.h"
 239 #include "AliGenPythiaEventHeader.h"
 240 #include "AliPythia.h"
 241 #include "AliRun.h"
 242
 243  ClassImp(AliGenPythia)
 244
 245 AliGenPythia::AliGenPythia()
 246                  :AliGenMC()
 247 {
 248 // Default Constructor
 249   fParticles = 0;
 250   fPythia    = 0;
 251   fDecayer   = new AliDecayerPythia();
 252   SetEventListRange();
 253   SetJetPhiRange();
 254   SetJetEtaRange();
 255   SetJetEtRange();
 256   SetGammaPhiRange();
 257   SetGammaEtaRange();
 258   SetPtKick();
 259 }
 260
 261 AliGenPythia::AliGenPythia(Int_t npart)
 262                  :AliGenMC(npart)
 263 {
 264 // default charm production at 5. 5 TeV
 265 // semimuonic decay
 266 // structure function GRVHO
 267 //
 268     fName = "Pythia";
 269     fTitle= "Particle Generator using PYTHIA";
 270     fXsection  = 0.;
 271     fNucA1=0;
 272     fNucA2=0;
 273     SetProcess();
 274     SetStrucFunc();
 275     SetForceDecay();
 276     SetPtHard();
 277     SetYHard();
 278     SetEnergyCMS();
 279     fDecayer = new AliDecayerPythia();
 280     // Set random number generator
 281     sRandom=fRandom;
 282     fFlavorSelect   = 0;
 283     // Produced particles
 284     fParticles = new TClonesArray("TParticle",1000);
 285     fEventVertex.Set(3);
 286     SetEventListRange();
 287     SetJetPhiRange();
 288     SetJetEtaRange();
 289     SetJetEtRange();
 290     SetGammaPhiRange();
 291     SetGammaEtaRange();
 292     SetJetReconstructionMode();
 293     SetPtKick();
 294     // Options determining what to keep in the stack (Heavy flavour generation)
 295     fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
 296     fFeedDownOpt = kTRUE;             // allow feed down from higher family
 297     // Fragmentation on/off
 298     fFragmentation = kTRUE;
 299     // Default counting mode
 300     fCountMode = kCountAll;
 301 }
 302
 303 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
 304 {
 305 // copy constructor
 306     Pythia.Copy(*this);
 307 }
 308
 309 AliGenPythia::~AliGenPythia()
 310 {
 311 // Destructor
 312 }
 313
 314 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
 315 {
 316   // Set a range of event numbers, for which a table
 317   // of generated particle will be printed
 318   fDebugEventFirst = eventFirst;
 319   fDebugEventLast  = eventLast;
 320   if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
 321 }
 322
 323 void AliGenPythia::Init()
 324 {
 325 // Initialisation
 326
 327     SetMC(AliPythia::Instance());
 328     fPythia=(AliPythia*) fgMCEvGen;
 329 //
 330     fParentWeight=1./Float_t(fNpart);
 331 //
 332 //  Forward Paramters to the AliPythia object
 333     fDecayer->SetForceDecay(fForceDecay);
 334     fDecayer->Init();
 335
 336
 337     fPythia->SetCKIN(3,fPtHardMin);
 338     fPythia->SetCKIN(4,fPtHardMax);
 339     fPythia->SetCKIN(7,fYHardMin);
 340     fPythia->SetCKIN(8,fYHardMax);
 341
 342     if (fNucA1 > 0 && fNucA2 > 0) fPythia->SetNuclei(fNucA1, fNucA2);
 343     // Fragmentation?
 344     if (fFragmentation) {
 345       fPythia->SetMSTP(111,1);
 346     } else {
 347       fPythia->SetMSTP(111,0);
 348     }
 349
 350
 351 //  initial state radiation
 352     fPythia->SetMSTP(61,fGinit);
 353 //  final state radiation
 354     fPythia->SetMSTP(71,fGfinal);
 355 //  pt - kick
 356     if (fPtKick > 0.) {
 357         fPythia->SetMSTP(91,1);
 358         fPythia->SetPARP(91,fPtKick);
 359     } else {
 360         fPythia->SetMSTP(91,0);
 361     }
 362
 363  //   fPythia->SetMSTJ(1,2);
 364  //
 365     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
 366
 367 //  Parent and Children Selection
 368     switch (fProcess)
 369     {
 370     case kPyCharm:
 371     case kPyCharmUnforced:
 372     case kPyCharmPbPbMNR:
 373     case kPyCharmppMNR:
 374     case kPyCharmpPbMNR:
 375         fParentSelect[0] =   411;
 376         fParentSelect[1] =   421;
 377         fParentSelect[2] =   431;
 378         fParentSelect[3] =  4122;
 379         fFlavorSelect    =  4;
 380         break;
 381     case kPyD0PbPbMNR:
 382     case kPyD0pPbMNR:
 383     case kPyD0ppMNR:
 384         fParentSelect[0] =   421;
 385         fFlavorSelect    =   4;
 386         break;
 387     case kPyBeauty:
 388     case kPyBeautyPbPbMNR:
 389     case kPyBeautypPbMNR:
 390     case kPyBeautyppMNR:
 391         fParentSelect[0]=  511;
 392         fParentSelect[1]=  521;
 393         fParentSelect[2]=  531;
 394         fParentSelect[3]= 5122;
 395         fParentSelect[4]= 5132;
 396         fParentSelect[5]= 5232;
 397         fParentSelect[6]= 5332;
 398         fFlavorSelect   = 5;
 399         break;
 400     case kPyBeautyUnforced:
 401         fParentSelect[0] =  511;
 402         fParentSelect[1] =  521;
 403         fParentSelect[2] =  531;
 404         fParentSelect[3] = 5122;
 405         fParentSelect[4] = 5132;
 406         fParentSelect[5] = 5232;
 407         fParentSelect[6] = 5332;
 408         fFlavorSelect    = 5;
 409         break;
 410     case kPyJpsiChi:
 411     case kPyJpsi:
 412         fParentSelect[0] = 443;
 413         break;
 414     case kPyMb:
 415     case kPyMbNonDiffr:
 416     case kPyJets:
 417     case kPyDirectGamma:
 418         break;
 419     }
 420 //
 421 //  This counts the total number of calls to Pyevnt() per run.
 422     fTrialsRun = 0;
 423     fQ         = 0.;
 424     fX1        = 0.;
 425     fX2        = 0.;
 426     fNev       = 0 ;
 427 //
 428     AliGenMC::Init();
 429 }
 430
 431 void AliGenPythia::Generate()
 432 {
 433 // Generate one event
 434     fDecayer->ForceDecay();
 435
 436     Float_t polar[3]   =   {0,0,0};
 437     Float_t origin[3]  =   {0,0,0};
 438     Float_t p[3];
 439 //  converts from mm/c to s
 440     const Float_t kconv=0.001/2.999792458e8;
 441 //
 442     Int_t nt=0;
 443     Int_t jev=0;
 444     Int_t j, kf;
 445     fTrials=0;
 446
 447     //  Set collision vertex position
 448     if(fVertexSmear==kPerEvent) {
 449         fPythia->SetMSTP(151,1);
 450         for (j=0;j<3;j++) {
 451             fPythia->SetPARP(151+j, fOsigma[j]*10.);
 452         }
 453     } else if (fVertexSmear==kPerTrack) {
 454         fPythia->SetMSTP(151,0);
 455     }
 456 //  event loop
 457     while(1)
 458     {
 459         fPythia->Pyevnt();
 460         if (gAlice->GetEvNumber()>=fDebugEventFirst &&
 461             gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
 462         fTrials++;
 463
 464         fPythia->ImportParticles(fParticles,"All");
 465
 466 //
 467 //
 468 //
 469         Int_t i;
 470
 471         Int_t np = fParticles->GetEntriesFast();
 472         if (np == 0 ) continue;
 473 // Get event vertex and discard the event if the Z coord. is too big
 474         TParticle *iparticle = (TParticle *) fParticles->At(0);
 475         Float_t distz = iparticle->Vz()/10.;
 476         if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue;
 477 //
 478         fEventVertex[0] = iparticle->Vx()/10.+fOrigin.At(0);
 479         fEventVertex[1] = iparticle->Vy()/10.+fOrigin.At(1);
 480         fEventVertex[2] = iparticle->Vz()/10.+fOrigin.At(2);
 481 //
 482         Int_t* pParent   = new Int_t[np];
 483         Int_t* pSelected = new Int_t[np];
 484         Int_t* trackIt   = new Int_t[np];
 485         for (i=0; i< np; i++) {
 486             pParent[i]   = -1;
 487             pSelected[i] =  0;
 488             trackIt[i]   =  0;
 489         }
 490
 491         Int_t nc = 0;        // Total n. of selected particles
 492         Int_t nParents = 0;  // Selected parents
 493         Int_t nTkbles = 0;   // Trackable particles
 494         if (fProcess != kPyMb && fProcess != kPyJets &&
 495             fProcess != kPyDirectGamma &&
 496             fProcess != kPyMbNonDiffr) {
 497
 498             for (i = 0; i<np; i++) {
 499                 iparticle = (TParticle *) fParticles->At(i);
 500                 Int_t ks = iparticle->GetStatusCode();
 501                 kf = CheckPDGCode(iparticle->GetPdgCode());
 502 // No initial state partons
 503                 if (ks==21) continue;
 504 //
 505 // Heavy Flavor Selection
 506 //
 507                 // quark ?
 508                 kf = TMath::Abs(kf);
 509                 Int_t kfl = kf;
 510                 // meson ?
 511                 if  (kfl > 10) kfl/=100;
 512                 // baryon
 513                 if (kfl > 10) kfl/=10;
 514                 if (kfl > 10) kfl/=10;
 515
 516                 Int_t ipa = iparticle->GetFirstMother()-1;
 517                 Int_t kfMo = 0;
 518
 519                 if (ipa > -1) {
 520                     TParticle *  mother = (TParticle *) fParticles->At(ipa);
 521                     kfMo = TMath::Abs(mother->GetPdgCode());
 522                 }
 523                 // What to keep in Stack?
 524                 Bool_t flavorOK = kFALSE;
 525                 Bool_t selectOK = kFALSE;
 526                 if (fFeedDownOpt) {
 527                   if (kfl >= fFlavorSelect) flavorOK = kTRUE;
 528                 } else {
 529                   if (kfl > fFlavorSelect) {
 530                     nc = -1;
 531                     break;
 532                   }
 533                   if (kfl == fFlavorSelect) flavorOK = kTRUE;
 534                 }
 535                 switch (fStackFillOpt) {
 536                 case kFlavorSelection:
 537                   selectOK = kTRUE;
 538                   break;
 539                 case kParentSelection:
 540                   if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
 541                   break;
 542                 }
 543                 if (flavorOK && selectOK) {
 544 //
 545 // Heavy flavor hadron or quark
 546 //
 547 // Kinematic seletion on final state heavy flavor mesons
 548                     if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
 549                     {
 550                       continue;
 551                     }
 552                     pSelected[i] = 1;
 553                     if (ParentSelected(kf)) ++nParents; // Update parent count
 554 //                  printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
 555                 } else {
 556 // Kinematic seletion on decay products
 557                     if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
 558                         && !KinematicSelection(iparticle, 1))
 559                     {
 560                       continue;
 561                     }
 562 //
 563 // Decay products
 564 // Select if mother was selected and is not tracked
 565
 566                     if (pSelected[ipa] &&
 567                         !trackIt[ipa]  &&     // mother will be  tracked ?
 568                         kfMo !=  5 &&         // mother is b-quark, don't store fragments
 569                         kfMo !=  4 &&         // mother is c-quark, don't store fragments
 570                         kf   != 92)           // don't store string
 571                     {
 572 //
 573 // Semi-stable or de-selected: diselect decay products:
 574 //
 575 //
 576                         if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
 577                         {
 578                             Int_t ipF = iparticle->GetFirstDaughter();
 579                             Int_t ipL = iparticle->GetLastDaughter();
 580                             if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
 581                         }
 582 //                      printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
 583                         pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
 584                     }
 585                 }
 586                 if (pSelected[i] == -1) pSelected[i] = 0;
 587                 if (!pSelected[i]) continue;
 588                 // Count quarks only if you did not include fragmentation
 589                 if (fFragmentation && kf <= 10) continue;
 590                 nc++;
 591 // Decision on tracking
 592                 trackIt[i] = 0;
 593 //
 594 // Track final state particle
 595                 if (ks == 1) trackIt[i] = 1;
 596 // Track semi-stable particles
 597                 if ((ks ==1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
 598 // Track particles selected by process if undecayed.
 599                 if (fForceDecay == kNoDecay) {
 600                     if (ParentSelected(kf)) trackIt[i] = 1;
 601                 } else {
 602                     if (ParentSelected(kf)) trackIt[i] = 0;
 603                 }
 604                 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
 605 //
 606 //
 607
 608             } // particle selection loop
 609             if (nc > 0) {
 610                 for (i = 0; i<np; i++) {
 611                     if (!pSelected[i]) continue;
 612                     TParticle *  iparticle = (TParticle *) fParticles->At(i);
 613                     kf = CheckPDGCode(iparticle->GetPdgCode());
 614                     Int_t ks = iparticle->GetStatusCode();
 615                     p[0] = iparticle->Px();
 616                     p[1] = iparticle->Py();
 617                     p[2] = iparticle->Pz();
 618                     origin[0] = fOrigin[0]+iparticle->Vx()/10.;
 619                     origin[1] = fOrigin[1]+iparticle->Vy()/10.;
 620                     origin[2] = fOrigin[2]+iparticle->Vz()/10.;
 621                     Float_t tof   = kconv*iparticle->T();
 622                     Int_t ipa     = iparticle->GetFirstMother()-1;
 623                     Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
 624                     SetTrack(fTrackIt*trackIt[i] ,
 625                                      iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1., ks);
 626                     pParent[i] = nt;
 627                     KeepTrack(nt);
 628                 } //  SetTrack loop
 629             }
 630         } else {
 631             nc = GenerateMB();
 632         } // mb ?
 633
 634         if (pParent)   delete[] pParent;
 635         if (pSelected) delete[] pSelected;
 636         if (trackIt)   delete[] trackIt;
 637
 638         if (nc > 0) {
 639           switch (fCountMode) {
 640           case kCountAll:
 641             // printf(" Count all \n");
 642             jev += nc;
 643             break;
 644           case kCountParents:
 645             // printf(" Count parents \n");
 646             jev += nParents;
 647             break;
 648           case kCountTrackables:
 649             // printf(" Count trackable \n");
 650             jev += nTkbles;
 651             break;
 652           }
 653             if (jev >= fNpart || fNpart == -1) {
 654                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 655                 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
 656
 657                 fQ  += fPythia->GetVINT(51);
 658                 fX1 += fPythia->GetVINT(41);
 659                 fX2 += fPythia->GetVINT(42);
 660                 fTrialsRun += fTrials;
 661                 fNev++;
 662                 MakeHeader();
 663                 break;
 664             }
 665         }
 666     } // event loop
 667     SetHighWaterMark(nt);
 668 //  adjust weight due to kinematic selection
 669     AdjustWeights();
 670 //  get cross-section
 671     fXsection=fPythia->GetPARI(1);
 672 }
 673
 674 Int_t  AliGenPythia::GenerateMB()
 675 {
 676 //
 677 // Min Bias selection and other global selections
 678 //
 679     Int_t i, kf, nt, iparent;
 680     Int_t nc = 0;
 681     Float_t p[3];
 682     Float_t polar[3]   =   {0,0,0};
 683     Float_t origin[3]  =   {0,0,0};
 684 //  converts from mm/c to s
 685     const Float_t kconv=0.001/2.999792458e8;
 686
 687     Int_t np = fParticles->GetEntriesFast();
 688     Int_t* pParent = new Int_t[np];
 689     for (i=0; i< np; i++) pParent[i] = -1;
 690     if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
 691         TParticle* jet1 = (TParticle *) fParticles->At(6);
 692         TParticle* jet2 = (TParticle *) fParticles->At(7);
 693         if (!CheckTrigger(jet1, jet2)) return 0;
 694     }
 695
 696     for (i = 0; i<np; i++) {
 697         Int_t trackIt = 0;
 698         TParticle *  iparticle = (TParticle *) fParticles->At(i);
 699         kf = CheckPDGCode(iparticle->GetPdgCode());
 700         Int_t ks = iparticle->GetStatusCode();
 701         Int_t km = iparticle->GetFirstMother();
 702         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
 703             (ks != 1) ||
 704             (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
 705             nc++;
 706             if (ks == 1) trackIt = 1;
 707             Int_t ipa = iparticle->GetFirstMother()-1;
 708
 709             iparent = (ipa > -1) ? pParent[ipa] : -1;
 710
 711 //
 712 // store track information
 713             p[0] = iparticle->Px();
 714             p[1] = iparticle->Py();
 715             p[2] = iparticle->Pz();
 716             origin[0] = fOrigin[0]+iparticle->Vx()/10.;
 717             origin[1] = fOrigin[1]+iparticle->Vy()/10.;
 718             origin[2] = fOrigin[2]+iparticle->Vz()/10.;
 719             Float_t tof=kconv*iparticle->T();
 720             SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
 721                      tof, kPPrimary, nt, 1., ks);
 722             KeepTrack(nt);
 723             pParent[i] = nt;
 724         } // select particle
 725     } // particle loop
 726
 727     if (pParent) delete[] pParent;
 728
 729     printf("\n I've put %i particles on the stack \n",nc);
 730     return nc;
 731 }
 732
 733
 734 void AliGenPythia::FinishRun()
 735 {
 736 // Print x-section summary
 737     fPythia->Pystat(1);
 738     fQ  /= fNev;
 739     fX1 /= fNev;
 740     fX2 /= fNev;
 741     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
 742     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
 743
 744
 745 }
 746
 747 void AliGenPythia::AdjustWeights()
 748 {
 749 // Adjust the weights after generation of all events
 750 //
 751     TParticle *part;
 752     Int_t ntrack=gAlice->GetNtrack();
 753     for (Int_t i=0; i<ntrack; i++) {
 754         part= gAlice->Particle(i);
 755         part->SetWeight(part->GetWeight()*fKineBias);
 756     }
 757 }
 758
 759 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
 760 {
 761 // Treat protons as inside nuclei with mass numbers a1 and a2
 762     fNucA1 = a1;
 763     fNucA2 = a2;
 764 }
 765
 766
 767 void AliGenPythia::MakeHeader()
 768 {
 769 // Builds the event header, to be called after each event
 770     AliGenEventHeader* header = new AliGenPythiaEventHeader("Pythia");
 771 //
 772 // Event type
 773     ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1));
 774 //
 775 // Number of trials
 776     ((AliGenPythiaEventHeader*) header)->SetTrials(fTrials);
 777 //
 778 // Event Vertex
 779     header->SetPrimaryVertex(fEventVertex);
 780 //
 781 // Jets that have triggered
 782     if (fProcess == kPyJets)
 783     {
 784         Int_t ntrig, njet;
 785         Float_t jets[4][10];
 786         GetJets(njet, ntrig, jets);
 787
 788         for (Int_t i = 0; i < ntrig; i++) {
 789             ((AliGenPythiaEventHeader*) header)->AddJet(jets[0][i], jets[1][i], jets[2][i],
 790                                                         jets[3][i]);
 791         }
 792     }
 793     gAlice->SetGenEventHeader(header);
 794 }
 795
 796
 797 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
 798 {
 799 // Check the kinematic trigger condition
 800 //
 801     Double_t eta[2];
 802     eta[0] = jet1->Eta();
 803     eta[1] = jet2->Eta();
 804     Double_t phi[2];
 805     phi[0] = jet1->Phi();
 806     phi[1] = jet2->Phi();
 807     Int_t    pdg[2];
 808     pdg[0] = jet1->GetPdgCode();
 809     pdg[1] = jet2->GetPdgCode();
 810     Bool_t   triggered = kFALSE;
 811
 812     if (fProcess == kPyJets) {
 813         Int_t njets = 0;
 814         Int_t ntrig = 0;
 815         Float_t jets[4][10];
 816 //
 817 // Use Pythia clustering on parton level to determine jet axis
 818 //
 819         GetJets(njets, ntrig, jets);
 820
 821         if (ntrig) triggered = kTRUE;
 822 //
 823     } else {
 824         Int_t ij = 0;
 825         Int_t ig = 1;
 826         if (pdg[0] == kGamma) {
 827             ij = 1;
 828             ig = 0;
 829         }
 830         //Check eta range first...
 831         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
 832             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
 833         {
 834             //Eta is okay, now check phi range
 835             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
 836                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
 837             {
 838                 triggered = kTRUE;
 839             }
 840         }
 841     }
 842     return triggered;
 843 }
 844
 845 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
 846 {
 847 // Assignment operator
 848     return *this;
 849 }
 850
 851 void  AliGenPythia::LoadEvent()
 852 {
 853 //
 854 // Load event into Pythia Common Block
 855 //
 856
 857
 858     Int_t npart = (Int_t) (gAlice->TreeK())->GetEntries();
 859    (fPythia->GetPyjets())->N = npart;
 860
 861     for (Int_t part = 0; part < npart; part++) {
 862         TParticle *MPart = gAlice->Particle(part);
 863         Int_t kf     = MPart->GetPdgCode();
 864         Int_t ks     = MPart->GetStatusCode();
 865         Float_t px = MPart->Px();
 866         Float_t py = MPart->Py();
 867         Float_t pz = MPart->Pz();
 868         Float_t e  = MPart->Energy();
 869         Float_t p  = TMath::Sqrt(px * px + py * py + pz * pz);
 870         Float_t m  = TMath::Sqrt(e * e - p * p);
 871
 872
 873         (fPythia->GetPyjets())->P[0][part] = px;
 874         (fPythia->GetPyjets())->P[1][part] = py;
 875         (fPythia->GetPyjets())->P[2][part] = pz;
 876         (fPythia->GetPyjets())->P[3][part] = e;
 877         (fPythia->GetPyjets())->P[4][part] = m;
 878
 879         (fPythia->GetPyjets())->K[1][part] = kf;
 880         (fPythia->GetPyjets())->K[0][part] = ks;
 881     }
 882 }
 883
 884 void AliGenPythia::RecJetsUA1(Float_t eCellMin, Float_t eCellSeed, Float_t eMin, Float_t rMax,
 885                               Int_t& njets, Float_t jets [4][50])
 886 {
 887 //
 888 //  Calls the Pythia jet finding algorithm to find jets in the current event
 889 //
 890 //
 891 //  Configure detector (EMCAL like)
 892 //
 893     fPythia->SetPARU(51,2.);
 894     fPythia->SetMSTU(51,Int_t(96 * 2./0.7));
 895     fPythia->SetMSTU(52,3 * 144);
 896 //
 897 //  Configure Jet Finder
 898 //
 899     fPythia->SetPARU(58, eCellMin);
 900     fPythia->SetPARU(52, eCellSeed);
 901     fPythia->SetPARU(53, eMin);
 902     fPythia->SetPARU(54, rMax);
 903     fPythia->SetMSTU(54, 2);
 904 //
 905 //  Save jets
 906     Int_t n     = fPythia->GetN();
 907
 908 //
 909 //  Run Jet Finder
 910     fPythia->Pycell(njets);
 911     Int_t i;
 912     for (i = 0; i < njets; i++) {
 913         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
 914         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
 915         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
 916         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
 917
 918         jets[0][i] = px;
 919         jets[1][i] = py;
 920         jets[2][i] = pz;
 921         jets[3][i] = e;
 922     }
 923 }
 924
 925
 926
 927 void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
 928 {
 929 //
 930 //  Calls the Pythia clustering algorithm to find jets in the current event
 931 //
 932     Int_t n     = fPythia->GetN();
 933     nJets       = 0;
 934     nJetsTrig   = 0;
 935     if (fJetReconstruction == kCluster) {
 936 //
 937 //  Configure cluster algorithm
 938 //
 939         fPythia->SetPARU(43, 2.);
 940         fPythia->SetMSTU(41, 1);
 941 //
 942 //  Call cluster algorithm
 943 //
 944         fPythia->Pyclus(nJets);
 945 //
 946 //  Loading jets from common block
 947 //
 948     } else {
 949 //
 950 //  Configure detector (EMCAL like)
 951 //
 952         fPythia->SetPARU(51,2.);
 953         fPythia->SetMSTU(51,Int_t(96 * 2./0.7));
 954         fPythia->SetMSTU(52,3 * 144);
 955 //
 956 //  Configure Jet Finder
 957 //
 958         fPythia->SetPARU(58,  0.0);
 959         fPythia->SetPARU(52,  4.0);
 960         fPythia->SetPARU(53, 10.0);
 961         fPythia->SetPARU(54,  1.0);
 962         fPythia->SetMSTU(54,  2);
 963 //
 964 //  Run Jet Finder
 965         fPythia->Pycell(nJets);
 966     }
 967
 968     Int_t i;
 969     for (i = 0; i < nJets; i++) {
 970         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
 971         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
 972         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
 973         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
 974         Float_t pt    = TMath::Sqrt(px * px + py * py);
 975         Float_t phi   = TMath::ATan2(py,px);
 976         Float_t theta = TMath::ATan2(pt,pz);
 977         Float_t et    = e * TMath::Sin(theta);
 978         Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
 979
 980         if (
 981             eta > fEtaMinJet && eta < fEtaMaxJet &&
 982             phi > fPhiMinJet && eta < fPhiMaxJet &&
 983             et  > fEtMinJet  && et  < fEtMaxJet
 984             )
 985         {
 986             jets[0][nJetsTrig] = px;
 987             jets[1][nJetsTrig] = py;
 988             jets[2][nJetsTrig] = pz;
 989             jets[3][nJetsTrig] = e;
 990             nJetsTrig++;
 991
 992         } else {
 993 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
 994         }
 995     }
 996 }
 997
 998
 999 #ifdef never
1000 void AliGenPythia::Streamer(TBuffer &R__b)
1001 {
1002    // Stream an object of class AliGenPythia.
1003
1004    if (R__b.IsReading()) {
1005       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1006       AliGenerator::Streamer(R__b);
1007       R__b >> (Int_t&)fProcess;
1008       R__b >> (Int_t&)fStrucFunc;
1009       R__b >> (Int_t&)fForceDecay;
1010       R__b >> fEnergyCMS;
1011       R__b >> fKineBias;
1012       R__b >> fTrials;
1013       fParentSelect.Streamer(R__b);
1014       fChildSelect.Streamer(R__b);
1015       R__b >> fXsection;
1016 //      (AliPythia::Instance())->Streamer(R__b);
1017       R__b >> fPtHardMin;
1018       R__b >> fPtHardMax;
1019 //      if (fDecayer) fDecayer->Streamer(R__b);
1020    } else {
1021       R__b.WriteVersion(AliGenPythia::IsA());
1022       AliGenerator::Streamer(R__b);
1023       R__b << (Int_t)fProcess;
1024       R__b << (Int_t)fStrucFunc;
1025       R__b << (Int_t)fForceDecay;
1026       R__b << fEnergyCMS;
1027       R__b << fKineBias;
1028       R__b << fTrials;
1029       fParentSelect.Streamer(R__b);
1030       fChildSelect.Streamer(R__b);
1031       R__b << fXsection;
1032 //      R__b << fPythia;
1033       R__b << fPtHardMin;
1034       R__b << fPtHardMax;
1035       //     fDecayer->Streamer(R__b);
1036    }
1037 }
1038 #endif
1039