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