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