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