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.41  2001/10/08 08:45:42  morsch
  19 Possibility of vertex cut added.
  20
  21 Revision 1.40  2001/09/25 11:30:23  morsch
  22 Pass event vertex to header.
  23
  24 Revision 1.39  2001/07/27 17:09:36  morsch
  25 Use local SetTrack, KeepTrack and SetHighWaterMark methods
  26 to delegate either to local stack or to stack owned by AliRun.
  27 (Piotr Skowronski, A.M.)
  28
  29 Revision 1.38  2001/07/13 10:58:54  morsch
  30 - Some coded moved to AliGenMC
  31 - Improved handling of secondary vertices.
  32
  33 Revision 1.37  2001/06/28 11:17:28  morsch
  34 SetEventListRange setter added. Events in specified range are listed for
  35 debugging. (Yuri Kharlov)
  36
  37 Revision 1.36  2001/03/30 07:05:49  morsch
  38 Final print-out in finish run.
  39 Write parton system for jet-production (preliminary solution).
  40
  41 Revision 1.35  2001/03/09 13:03:40  morsch
  42 Process_t and Struc_Func_t moved to AliPythia.h
  43
  44 Revision 1.34  2001/02/14 15:50:40  hristov
  45 The last particle in event marked using SetHighWaterMark
  46
  47 Revision 1.33  2001/01/30 09:23:12  hristov
  48 Streamers removed (R.Brun)
  49
  50 Revision 1.32  2001/01/26 19:55:51  hristov
  51 Major upgrade of AliRoot code
  52
  53 Revision 1.31  2001/01/17 10:54:31  hristov
  54 Better protection against FPE
  55
  56 Revision 1.30  2000/12/18 08:55:35  morsch
  57 Make AliPythia dependent generartors work with new scheme of random number generation
  58
  59 Revision 1.29  2000/12/04 11:22:03  morsch
  60 Init of sRandom as in 1.15
  61
  62 Revision 1.28  2000/12/02 11:41:39  morsch
  63 Use SetRandom() to initialize random number generator in constructor.
  64
  65 Revision 1.27  2000/11/30 20:29:02  morsch
  66 Initialise static variable sRandom in constructor: sRandom = fRandom;
  67
  68 Revision 1.26  2000/11/30 07:12:50  alibrary
  69 Introducing new Rndm and QA classes
  70
  71 Revision 1.25  2000/10/18 19:11:27  hristov
  72 Division by zero fixed
  73
  74 Revision 1.24  2000/09/18 10:41:35  morsch
  75 Add possibility to use nuclear structure functions from PDF library V8.
  76
  77 Revision 1.23  2000/09/14 14:05:40  morsch
  78 dito
  79
  80 Revision 1.22  2000/09/14 14:02:22  morsch
  81 - Correct conversion from mm to cm when passing particle vertex to MC.
  82 - Correct handling of fForceDecay == all.
  83
  84 Revision 1.21  2000/09/12 14:14:55  morsch
  85 Call fDecayer->ForceDecay() at the beginning of Generate().
  86
  87 Revision 1.20  2000/09/06 14:29:33  morsch
  88 Use AliPythia for event generation an AliDecayPythia for decays.
  89 Correct handling of "nodecay" option
  90
  91 Revision 1.19  2000/07/11 18:24:56  fca
  92 Coding convention corrections + few minor bug fixes
  93
  94 Revision 1.18  2000/06/30 12:40:34  morsch
  95 Pythia takes care of vertex smearing. Correct conversion from Pythia units (mm) to
  96 Geant units (cm).
  97
  98 Revision 1.17  2000/06/09 20:34:07  morsch
  99 All coding rule violations except RS3 corrected
 100
 101 Revision 1.16  2000/05/15 15:04:20  morsch
 102 The full event is written for fNtrack = -1
 103 Coding rule violations corrected.
 104
 105 Revision 1.15  2000/04/26 10:14:24  morsch
 106 Particles array has one entry more than pythia particle list. Upper bound of
 107 particle loop changed to np-1 (R. Guernane, AM)
 108
 109 Revision 1.14  2000/04/05 08:36:13  morsch
 110 Check status code of particles in Pythia event
 111 to avoid double counting as partonic state and final state particle.
 112
 113 Revision 1.13  1999/11/09 07:38:48  fca
 114 Changes for compatibility with version 2.23 of ROOT
 115
 116 Revision 1.12  1999/11/03 17:43:20  fca
 117 New version from G.Martinez & A.Morsch
 118
 119 Revision 1.11  1999/09/29 09:24:14  fca
 120 Introduction of the Copyright and cvs Log
 121 */
 122
 123 #include "AliGenPythia.h"
 124 #include "AliGenPythiaEventHeader.h"
 125 #include "AliDecayerPythia.h"
 126 #include "AliRun.h"
 127 #include "AliPythia.h"
 128 #include "AliPDG.h"
 129 #include <TParticle.h>
 130 #include <TSystem.h>
 131
 132  ClassImp(AliGenPythia)
 133
 134 AliGenPythia::AliGenPythia()
 135                  :AliGenMC()
 136 {
 137 // Default Constructor
 138   fDecayer = new AliDecayerPythia();
 139   SetEventListRange();
 140 }
 141
 142 AliGenPythia::AliGenPythia(Int_t npart)
 143                  :AliGenMC(npart)
 144 {
 145 // default charm production at 5. 5 TeV
 146 // semimuonic decay
 147 // structure function GRVHO
 148 //
 149     fXsection  = 0.;
 150     fNucA1=0;
 151     fNucA2=0;
 152     SetProcess();
 153     SetStrucFunc();
 154     SetForceDecay();
 155     SetPtHard();
 156     SetEnergyCMS();
 157     fDecayer = new AliDecayerPythia();
 158     // Set random number generator
 159     sRandom=fRandom;
 160     SetEventListRange();
 161     fFlavorSelect   = 0;
 162     // Produced particles
 163     fParticles = new TClonesArray("TParticle",1000);
 164     fEventVertex.Set(3);
 165 }
 166
 167 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
 168 {
 169 // copy constructor
 170 }
 171
 172 AliGenPythia::~AliGenPythia()
 173 {
 174 // Destructor
 175 }
 176
 177 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
 178 {
 179   // Set a range of event numbers, for which a table
 180   // of generated particle will be printed
 181   fDebugEventFirst = eventFirst;
 182   fDebugEventLast  = eventLast;
 183   if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
 184 }
 185
 186 void AliGenPythia::Init()
 187 {
 188 // Initialisation
 189   SetMC(AliPythia::Instance());
 190     fPythia=(AliPythia*) fgMCEvGen;
 191 //
 192     fParentWeight=1./Float_t(fNpart);
 193 //
 194 //  Forward Paramters to the AliPythia object
 195     fDecayer->SetForceDecay(fForceDecay);
 196     fDecayer->Init();
 197
 198
 199     fPythia->SetCKIN(3,fPtHardMin);
 200     fPythia->SetCKIN(4,fPtHardMax);
 201     if (fNucA1 > 0 && fNucA2 > 0) fPythia->SetNuclei(fNucA1, fNucA2);
 202     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
 203
 204     //    fPythia->Pylist(0);
 205     //    fPythia->Pystat(2);
 206 //  Parent and Children Selection
 207     switch (fProcess)
 208     {
 209     case kPyCharm:
 210         fParentSelect[0] =  411;
 211         fParentSelect[1] =  421;
 212         fParentSelect[2] =  431;
 213         fParentSelect[3] = 4122;
 214         fFlavorSelect    = 4;
 215         break;
 216     case kPyCharmUnforced:
 217         fParentSelect[0] =  411;
 218         fParentSelect[1] =  421;
 219         fParentSelect[2] =  431;
 220         fParentSelect[3]=  4122;
 221         fFlavorSelect    = 4;
 222         break;
 223     case kPyBeauty:
 224         fParentSelect[0]=  511;
 225         fParentSelect[1]=  521;
 226         fParentSelect[2]=  531;
 227         fParentSelect[3]= 5122;
 228         fParentSelect[4]= 5132;
 229         fParentSelect[5]= 5232;
 230         fParentSelect[6]= 5332;
 231         fFlavorSelect   = 5;
 232         break;
 233     case kPyBeautyUnforced:
 234         fParentSelect[0] =  511;
 235         fParentSelect[1] =  521;
 236         fParentSelect[2] =  531;
 237         fParentSelect[3] = 5122;
 238         fParentSelect[4] = 5132;
 239         fParentSelect[5] = 5232;
 240         fParentSelect[6] = 5332;
 241         fFlavorSelect    = 5;
 242         break;
 243     case kPyJpsiChi:
 244     case kPyJpsi:
 245         fParentSelect[0] = 443;
 246         break;
 247     case kPyMb:
 248     case kPyJets:
 249     case kPyDirectGamma:
 250         break;
 251     }
 252     AliGenMC::Init();
 253 }
 254
 255 void AliGenPythia::Generate()
 256 {
 257 // Generate one event
 258     fDecayer->ForceDecay();
 259
 260     Float_t polar[3]   =   {0,0,0};
 261     Float_t origin[3]  =   {0,0,0};
 262     Float_t p[3];
 263 //  converts from mm/c to s
 264     const Float_t kconv=0.001/2.999792458e8;
 265 //
 266     Int_t nt=0;
 267     Int_t jev=0;
 268     Int_t j, kf;
 269     fTrials=0;
 270
 271 //  Set collision vertex position
 272     if(fVertexSmear==kPerEvent) {
 273         fPythia->SetMSTP(151,1);
 274         for (j=0;j<3;j++) {
 275             fPythia->SetPARP(151+j, fOsigma[j]*10.);
 276         }
 277     } else if (fVertexSmear==kPerTrack) {
 278         fPythia->SetMSTP(151,0);
 279     }
 280 //  event loop
 281     while(1)
 282     {
 283         fPythia->Pyevnt();
 284         if (gAlice->GetEvNumber()>=fDebugEventFirst &&
 285             gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
 286         fTrials++;
 287         fPythia->ImportParticles(fParticles,"All");
 288
 289 //
 290 //
 291 //
 292         Int_t i;
 293
 294         Int_t np = fParticles->GetEntriesFast();
 295         if (np == 0 ) continue;
 296 // Get event vertex and discard the event if the Z coord. is too big
 297         TParticle *iparticle = (TParticle *) fParticles->At(0);
 298         Float_t distz = iparticle->Vz()/10.;
 299         if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue;
 300 //
 301         fEventVertex[0] = iparticle->Vx()/10.+fOrigin.At(0);
 302         fEventVertex[1] = iparticle->Vy()/10.+fOrigin.At(1);
 303         fEventVertex[2] = iparticle->Vz()/10.+fOrigin.At(2);
 304 //
 305         Int_t* pParent   = new Int_t[np];
 306         Int_t* pSelected = new Int_t[np];
 307         Int_t* trackIt   = new Int_t[np];
 308         for (i=0; i< np-1; i++) {
 309             pParent[i]   = -1;
 310             pSelected[i] =  0;
 311         }
 312         printf("\n **************************************************%d\n",np);
 313         Int_t nc = 0;
 314         if (fProcess != kPyMb && fProcess != kPyJets && fProcess != kPyDirectGamma) {
 315             for (i = 0; i<np-1; i++) {
 316                 iparticle = (TParticle *) fParticles->At(i);
 317                 Int_t ks = iparticle->GetStatusCode();
 318                 kf = CheckPDGCode(iparticle->GetPdgCode());
 319 // No initial state partons
 320                 if (ks==21) continue;
 321 //
 322 // Heavy Flavor Selection
 323 //
 324                 // quark ?
 325                 kf = TMath::Abs(kf);
 326                 Int_t kfl = kf;
 327                 // meson ?
 328                 if  (kfl > 10) kfl/=100;
 329                 // baryon
 330                 if (kfl > 10) kfl/=10;
 331                 if (kfl > 10) kfl/=10;
 332
 333                 Int_t ipa = iparticle->GetFirstMother()-1;
 334                 Int_t kfMo = 0;
 335
 336                 if (ipa > -1) {
 337                     TParticle *  mother = (TParticle *) fParticles->At(ipa);
 338                     kfMo = TMath::Abs(mother->GetPdgCode());
 339                 }
 340 //              printf("\n particle (all)  %d %d %d", i, pSelected[i], kf);
 341                 if (kfl >= fFlavorSelect) {
 342 //
 343 // Heavy flavor hadron or quark
 344 //
 345 // Kinematic seletion on final state heavy flavor mesons
 346                     if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
 347                     {
 348                         nc = -1;
 349                         break;
 350                     }
 351                     pSelected[i] = 1;
 352 //                  printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
 353                 } else {
 354 // Kinematic seletion on decay products
 355                     if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
 356                         && !KinematicSelection(iparticle, 1))
 357                     {
 358                         nc = -1;
 359                         break;
 360                     }
 361 //
 362 // Decay products
 363 // Select if mother was selected and is not tracked
 364
 365                     if (pSelected[ipa] &&
 366                         !trackIt[ipa]  &&     // mother will be  tracked ?
 367                         kfMo !=  5 &&         // mother is b-quark, don't store fragments
 368                         kfMo !=  4 &&         // mother is c-quark, don't store fragments
 369                         kf   != 92)           // don't store string
 370                     {
 371 //
 372 // Semi-stable or de-selected: diselect decay products:
 373 //
 374 //
 375                         if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > 10e-15)
 376                         {
 377                             Int_t ipF = iparticle->GetFirstDaughter();
 378                             Int_t ipL = iparticle->GetLastDaughter();
 379                             if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
 380                         }
 381 //                      printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
 382                         pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
 383                     }
 384                 }
 385                 if (pSelected[i] == -1) pSelected[i] = 0;
 386                 if (!pSelected[i]) continue;
 387                 nc++;
 388 // Decision on tracking
 389                 trackIt[i] = 0;
 390 //
 391 // Track final state particle
 392                 if (ks == 1) trackIt[i] = 1;
 393 // Track semi-stable particles
 394                 if ((ks ==1) || (fDecayer->GetLifetime(kf)> 10e-15))  trackIt[i] = 1;
 395 // Track particles selected by process if undecayed.
 396                 if (fForceDecay == kNoDecay) {
 397                     if (ParentSelected(kf)) trackIt[i] = 1;
 398                 } else {
 399                     if (ParentSelected(kf)) trackIt[i] = 0;
 400                 }
 401 //
 402 //
 403
 404             } // particle selection loop
 405             if (nc > -1) {
 406                 for (i = 0; i<np-1; i++) {
 407                     if (!pSelected[i]) continue;
 408                     TParticle *  iparticle = (TParticle *) fParticles->At(i);
 409                     kf = CheckPDGCode(iparticle->GetPdgCode());
 410                     p[0] = iparticle->Px();
 411                     p[1] = iparticle->Py();
 412                     p[2] = iparticle->Pz();
 413                     origin[0] = fOrigin[0]+iparticle->Vx()/10.;
 414                     origin[1] = fOrigin[1]+iparticle->Vy()/10.;
 415                     origin[2] = fOrigin[2]+iparticle->Vz()/10.;
 416                     Float_t tof   = kconv*iparticle->T();
 417                     Int_t ipa     = iparticle->GetFirstMother()-1;
 418                     Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
 419                     SetTrack(fTrackIt*trackIt[i] ,
 420                                      iparent, kf, p, origin, polar, tof, kPPrimary, nt, 1.);
 421                     pParent[i] = nt;
 422                     KeepTrack(nt);
 423                 } //  SetTrack loop
 424             }
 425         } else {
 426             nc = GenerateMB();
 427         } // mb ?
 428
 429         if (nc > 0) {
 430             jev+=nc;
 431             if (jev >= fNpart || fNpart == -1) {
 432                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 433                 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
 434                 break;
 435             }
 436         }
 437     } // event loop
 438     SetHighWaterMark(nt);
 439 //  adjust weight due to kinematic selection
 440     AdjustWeights();
 441 //  get cross-section
 442     fXsection=fPythia->GetPARI(1);
 443 }
 444
 445 Int_t  AliGenPythia::GenerateMB()
 446 {
 447     Int_t i, kf, nt, iparent;
 448     Int_t nc = 0;
 449     Float_t p[3];
 450     Float_t polar[3]   =   {0,0,0};
 451     Float_t origin[3]  =   {0,0,0};
 452 //  converts from mm/c to s
 453     const Float_t kconv=0.001/2.999792458e8;
 454
 455     Int_t np = fParticles->GetEntriesFast();
 456     Int_t* pParent = new Int_t[np];
 457     for (i=0; i< np-1; i++) pParent[i] = -1;
 458
 459     for (i = 0; i<np-1; i++) {
 460         Int_t trackIt = 0;
 461         TParticle *  iparticle = (TParticle *) fParticles->At(i);
 462         kf = CheckPDGCode(iparticle->GetPdgCode());
 463         Int_t ks = iparticle->GetStatusCode();
 464         Int_t km = iparticle->GetFirstMother();
 465 //      printf("\n Particle: %d %d %d", i, kf, ks);
 466
 467         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
 468             (ks != 1) ||
 469             (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
 470             nc++;
 471             if (ks == 1) trackIt = 1;
 472             Int_t ipa = iparticle->GetFirstMother()-1;
 473
 474             iparent = (ipa > -1) ? pParent[ipa] : -1;
 475
 476 //
 477 // store track information
 478             p[0] = iparticle->Px();
 479             p[1] = iparticle->Py();
 480             p[2] = iparticle->Pz();
 481             origin[0] = fOrigin[0]+iparticle->Vx()/10.;
 482             origin[1] = fOrigin[1]+iparticle->Vy()/10.;
 483             origin[2] = fOrigin[2]+iparticle->Vz()/10.;
 484             Float_t tof=kconv*iparticle->T();
 485             SetTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
 486                          tof, kPPrimary, nt);
 487             KeepTrack(nt);
 488             pParent[i] = nt;
 489         } // select particle
 490     } // particle loop
 491
 492     if (pParent) delete[] pParent;
 493
 494     printf("\n I've put %i particles on the stack \n",nc);
 495     MakeHeader();
 496     return nc;
 497 }
 498
 499
 500 void AliGenPythia::FinishRun()
 501 {
 502 // Print x-section summary
 503     fPythia->Pystat(1);
 504 }
 505
 506 void AliGenPythia::AdjustWeights()
 507 {
 508 // Adjust the weights after generation of all events
 509 //
 510     TParticle *part;
 511     Int_t ntrack=gAlice->GetNtrack();
 512     for (Int_t i=0; i<ntrack; i++) {
 513         part= gAlice->Particle(i);
 514         part->SetWeight(part->GetWeight()*fKineBias);
 515     }
 516 }
 517
 518 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
 519 {
 520 // Treat protons as inside nuclei with mass numbers a1 and a2
 521     fNucA1 = a1;
 522     fNucA2 = a2;
 523 }
 524
 525
 526 void AliGenPythia::MakeHeader()
 527 {
 528 // Builds the event header, to be called after each event
 529     AliGenEventHeader* header = new AliGenPythiaEventHeader("Pythia");
 530     ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1));
 531     header->SetPrimaryVertex(fEventVertex);
 532     gAlice->SetGenEventHeader(header);
 533 }
 534
 535
 536 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
 537 {
 538 // Assignment operator
 539     return *this;
 540 }
 541
 542
 543
 544 #ifdef never
 545 void AliGenPythia::Streamer(TBuffer &R__b)
 546 {
 547    // Stream an object of class AliGenPythia.
 548
 549    if (R__b.IsReading()) {
 550       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
 551       AliGenerator::Streamer(R__b);
 552       R__b >> (Int_t&)fProcess;
 553       R__b >> (Int_t&)fStrucFunc;
 554       R__b >> (Int_t&)fForceDecay;
 555       R__b >> fEnergyCMS;
 556       R__b >> fKineBias;
 557       R__b >> fTrials;
 558       fParentSelect.Streamer(R__b);
 559       fChildSelect.Streamer(R__b);
 560       R__b >> fXsection;
 561 //      (AliPythia::Instance())->Streamer(R__b);
 562       R__b >> fPtHardMin;
 563       R__b >> fPtHardMax;
 564 //      if (fDecayer) fDecayer->Streamer(R__b);
 565    } else {
 566       R__b.WriteVersion(AliGenPythia::IsA());
 567       AliGenerator::Streamer(R__b);
 568       R__b << (Int_t)fProcess;
 569       R__b << (Int_t)fStrucFunc;
 570       R__b << (Int_t)fForceDecay;
 571       R__b << fEnergyCMS;
 572       R__b << fKineBias;
 573       R__b << fTrials;
 574       fParentSelect.Streamer(R__b);
 575       fChildSelect.Streamer(R__b);
 576       R__b << fXsection;
 577 //      R__b << fPythia;
 578       R__b << fPtHardMin;
 579       R__b << fPtHardMax;
 580       //     fDecayer->Streamer(R__b);
 581    }
 582 }
 583 #endif
 584