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