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