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