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