PYTHIA6/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 /* $Id$ */
  17
  18 //
  19 // Generator using the TPythia interface (via AliPythia)
  20 // to generate pp collisions.
  21 // Using SetNuclei() also nuclear modifications to the structure functions
  22 // can be taken into account. This makes, of course, only sense for the
  23 // generation of the products of hard processes (heavy flavor, jets ...)
  24 //
  25 // andreas.morsch@cern.ch
  26 //
  27
  28 #include <TDatabasePDG.h>
  29 #include <TParticle.h>
  30 #include <TPDGCode.h>
  31 #include <TSystem.h>
  32 #include <TTree.h>
  33
  34 #include "AliConst.h"
  35 #include "AliDecayerPythia.h"
  36 #include "AliGenPythia.h"
  37 #include "AliHeader.h"
  38 #include "AliGenPythiaEventHeader.h"
  39 #include "AliPythia.h"
  40 #include "AliPythiaRndm.h"
  41 #include "AliRun.h"
  42 #include "AliStack.h"
  43 #include "AliRunLoader.h"
  44 #include "AliMC.h"
  45
  46 ClassImp(AliGenPythia)
  47
  48 AliGenPythia::AliGenPythia()
  49                  :AliGenMC()
  50 {
  51 // Default Constructor
  52   fParticles = 0;
  53   fPythia    = 0;
  54   fHeader = 0;
  55   fDecayer   = new AliDecayerPythia();
  56   SetEventListRange();
  57   SetJetPhiRange();
  58   SetJetEtaRange();
  59   SetJetEtRange();
  60   SetGammaPhiRange();
  61   SetGammaEtaRange();
  62   SetPtKick();
  63   SetQuench();
  64   SetHadronisation();
  65   fSetNuclei = kFALSE;
  66   if (!AliPythiaRndm::GetPythiaRandom())
  67     AliPythiaRndm::SetPythiaRandom(GetRandom());
  68 }
  69
  70 AliGenPythia::AliGenPythia(Int_t npart)
  71                  :AliGenMC(npart)
  72 {
  73 // default charm production at 5. 5 TeV
  74 // semimuonic decay
  75 // structure function GRVHO
  76 //
  77     fName = "Pythia";
  78     fTitle= "Particle Generator using PYTHIA";
  79     fXsection  = 0.;
  80     SetProcess();
  81     SetStrucFunc();
  82     SetForceDecay();
  83     SetPtHard();
  84     SetYHard();
  85     SetEnergyCMS();
  86     fDecayer = new AliDecayerPythia();
  87     // Set random number generator
  88     if (!AliPythiaRndm::GetPythiaRandom())
  89       AliPythiaRndm::SetPythiaRandom(GetRandom());
  90     fFlavorSelect   = 0;
  91     // Produced particles
  92     fParticles = new TClonesArray("TParticle",1000);
  93     fHeader = 0;
  94     SetEventListRange();
  95     SetJetPhiRange();
  96     SetJetEtaRange();
  97     SetJetEtRange();
  98     SetGammaPhiRange();
  99     SetGammaEtaRange();
 100     SetJetReconstructionMode();
 101     SetQuench();
 102     SetHadronisation();
 103     SetPtKick();
 104     // Options determining what to keep in the stack (Heavy flavour generation)
 105     fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
 106     fFeedDownOpt = kTRUE;             // allow feed down from higher family
 107     // Fragmentation on/off
 108     fFragmentation = kTRUE;
 109     // Default counting mode
 110     fCountMode = kCountAll;
 111     // Pycel
 112     SetPycellParameters();
 113     fSetNuclei = kFALSE;
 114 }
 115
 116 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
 117     :AliGenMC(Pythia)
 118 {
 119 // copy constructor
 120     Pythia.Copy(*this);
 121 }
 122
 123 AliGenPythia::~AliGenPythia()
 124 {
 125 // Destructor
 126 }
 127
 128 void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
 129                                        Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
 130 {
 131 // Set pycell parameters
 132     fPycellEtaMax    =  etamax;
 133     fPycellNEta      =  neta;
 134     fPycellNPhi      =  nphi;
 135     fPycellThreshold =  thresh;
 136     fPycellEtSeed    =  etseed;
 137     fPycellMinEtJet  =  minet;
 138     fPycellMaxRadius =  r;
 139 }
 140
 141
 142
 143 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
 144 {
 145   // Set a range of event numbers, for which a table
 146   // of generated particle will be printed
 147   fDebugEventFirst = eventFirst;
 148   fDebugEventLast  = eventLast;
 149   if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
 150 }
 151
 152 void AliGenPythia::Init()
 153 {
 154 // Initialisation
 155
 156     SetMC(AliPythia::Instance());
 157     fPythia=(AliPythia*) fMCEvGen;
 158
 159 //
 160     fParentWeight=1./Float_t(fNpart);
 161 //
 162 //  Forward Paramters to the AliPythia object
 163     fDecayer->SetForceDecay(fForceDecay);
 164     fDecayer->Init();
 165
 166
 167     fPythia->SetCKIN(3,fPtHardMin);
 168     fPythia->SetCKIN(4,fPtHardMax);
 169     fPythia->SetCKIN(7,fYHardMin);
 170     fPythia->SetCKIN(8,fYHardMax);
 171
 172     if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget);
 173     // Fragmentation?
 174     if (fFragmentation) {
 175       fPythia->SetMSTP(111,1);
 176     } else {
 177       fPythia->SetMSTP(111,0);
 178     }
 179
 180
 181 //  initial state radiation
 182     fPythia->SetMSTP(61,fGinit);
 183 //  final state radiation
 184     fPythia->SetMSTP(71,fGfinal);
 185 //  pt - kick
 186     if (fPtKick > 0.) {
 187         fPythia->SetMSTP(91,1);
 188         fPythia->SetPARP(91,fPtKick);
 189     } else {
 190         fPythia->SetMSTP(91,0);
 191     }
 192
 193
 194     if (fReadFromFile) {
 195         fRL  =  AliRunLoader::Open(fFileName, "Partons");
 196         fRL->LoadKinematics();
 197         fRL->LoadHeader();
 198     } else {
 199         fRL = 0x0;
 200     }
 201
 202
 203  //
 204     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
 205
 206 //  Parent and Children Selection
 207     switch (fProcess)
 208     {
 209     case kPyCharm:
 210     case kPyCharmUnforced:
 211     case kPyCharmPbPbMNR:
 212     case kPyCharmppMNR:
 213     case kPyCharmpPbMNR:
 214         fParentSelect[0] =   411;
 215         fParentSelect[1] =   421;
 216         fParentSelect[2] =   431;
 217         fParentSelect[3] =  4122;
 218         fFlavorSelect    =  4;
 219         break;
 220     case kPyD0PbPbMNR:
 221     case kPyD0pPbMNR:
 222     case kPyD0ppMNR:
 223         fParentSelect[0] =   421;
 224         fFlavorSelect    =   4;
 225         break;
 226     case kPyDPlusPbPbMNR:
 227     case kPyDPluspPbMNR:
 228     case kPyDPlusppMNR:
 229         fParentSelect[0] =   411;
 230         fFlavorSelect    =   4;
 231         break;
 232     case kPyBeauty:
 233     case kPyBeautyPbPbMNR:
 234     case kPyBeautypPbMNR:
 235     case kPyBeautyppMNR:
 236         fParentSelect[0]=  511;
 237         fParentSelect[1]=  521;
 238         fParentSelect[2]=  531;
 239         fParentSelect[3]= 5122;
 240         fParentSelect[4]= 5132;
 241         fParentSelect[5]= 5232;
 242         fParentSelect[6]= 5332;
 243         fFlavorSelect   = 5;
 244         break;
 245     case kPyBeautyUnforced:
 246         fParentSelect[0] =  511;
 247         fParentSelect[1] =  521;
 248         fParentSelect[2] =  531;
 249         fParentSelect[3] = 5122;
 250         fParentSelect[4] = 5132;
 251         fParentSelect[5] = 5232;
 252         fParentSelect[6] = 5332;
 253         fFlavorSelect    = 5;
 254         break;
 255     case kPyJpsiChi:
 256     case kPyJpsi:
 257         fParentSelect[0] = 443;
 258         break;
 259     case kPyMb:
 260     case kPyMbNonDiffr:
 261     case kPyJets:
 262     case kPyDirectGamma:
 263         break;
 264     }
 265 //
 266 //
 267 //  JetFinder for Trigger
 268 //
 269 //  Configure detector (EMCAL like)
 270 //
 271         fPythia->SetPARU(51, fPycellEtaMax);
 272         fPythia->SetMSTU(51, fPycellNEta);
 273         fPythia->SetMSTU(52, fPycellNPhi);
 274 //
 275 //  Configure Jet Finder
 276 //
 277         fPythia->SetPARU(58,  fPycellThreshold);
 278         fPythia->SetPARU(52,  fPycellEtSeed);
 279         fPythia->SetPARU(53,  fPycellMinEtJet);
 280         fPythia->SetPARU(54,  fPycellMaxRadius);
 281         fPythia->SetMSTU(54,  2);
 282 //
 283 //  This counts the total number of calls to Pyevnt() per run.
 284     fTrialsRun = 0;
 285     fQ         = 0.;
 286     fX1        = 0.;
 287     fX2        = 0.;
 288     fNev       = 0 ;
 289 //
 290 //
 291 //
 292     AliGenMC::Init();
 293 //
 294 //
 295 //
 296     if (fSetNuclei) {
 297         fDyBoost = 0;
 298         Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
 299     }
 300
 301     if (fQuench) {
 302         fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
 303     }
 304
 305 }
 306
 307 void AliGenPythia::Generate()
 308 {
 309 // Generate one event
 310
 311     fDecayer->ForceDecay();
 312
 313     Float_t polar[3]   =   {0,0,0};
 314     Float_t origin[3]  =   {0,0,0};
 315     Float_t p[4];
 316 //  converts from mm/c to s
 317     const Float_t kconv=0.001/2.999792458e8;
 318 //
 319     Int_t nt=0;
 320     Int_t jev=0;
 321     Int_t j, kf;
 322     fTrials=0;
 323
 324
 325     //  Set collision vertex position
 326     if (fVertexSmear == kPerEvent) Vertex();
 327
 328 //  event loop
 329     while(1)
 330     {
 331 //
 332 // Produce event
 333 //
 334 //
 335 // Switch hadronisation off
 336 //
 337         fPythia->SetMSTJ(1, 0);
 338 //
 339 // Either produce new event or read partons from file
 340 //
 341         if (!fReadFromFile) {
 342             fPythia->Pyevnt();
 343             fNpartons = fPythia->GetN();
 344         } else {
 345             printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
 346             fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
 347             fPythia->SetN(0);
 348             LoadEvent(fRL->Stack(), 0 , 1);
 349             fPythia->Pyedit(21);
 350         }
 351
 352 //
 353 //  Run quenching routine
 354 //
 355         if (fQuench == 1) {
 356             fPythia->Quench();
 357         } else if (fQuench == 2){
 358             fPythia->Pyquen(208., 0, 0.);
 359         }
 360 //
 361 // Switch hadronisation on
 362 //
 363         fPythia->SetMSTJ(1, 1);
 364 //
 365 // .. and perform hadronisation
 366 //      printf("Calling hadronisation %d\n", fPythia->GetN());
 367         fPythia->Pyexec();
 368
 369         if (gAlice) {
 370             if (gAlice->GetEvNumber()>=fDebugEventFirst &&
 371                 gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
 372         }
 373
 374         fTrials++;
 375         fPythia->ImportParticles(fParticles,"All");
 376         Boost();
 377 //
 378 //
 379 //
 380         Int_t i;
 381
 382
 383         Int_t np = fParticles->GetEntriesFast();
 384
 385         if (np == 0 ) continue;
 386 //
 387
 388 //
 389         Int_t* pParent   = new Int_t[np];
 390         Int_t* pSelected = new Int_t[np];
 391         Int_t* trackIt   = new Int_t[np];
 392         for (i = 0; i < np; i++) {
 393             pParent[i]   = -1;
 394             pSelected[i] =  0;
 395             trackIt[i]   =  0;
 396         }
 397
 398         Int_t nc = 0;        // Total n. of selected particles
 399         Int_t nParents = 0;  // Selected parents
 400         Int_t nTkbles = 0;   // Trackable particles
 401         if (fProcess != kPyMb && fProcess != kPyJets &&
 402             fProcess != kPyDirectGamma &&
 403             fProcess != kPyMbNonDiffr) {
 404
 405             for (i = 0; i < np; i++) {
 406                 TParticle* iparticle = (TParticle *) fParticles->At(i);
 407                 Int_t ks = iparticle->GetStatusCode();
 408                 kf = CheckPDGCode(iparticle->GetPdgCode());
 409 // No initial state partons
 410                 if (ks==21) continue;
 411 //
 412 // Heavy Flavor Selection
 413 //
 414                 // quark ?
 415                 kf = TMath::Abs(kf);
 416                 Int_t kfl = kf;
 417                 // meson ?
 418                 if  (kfl > 10) kfl/=100;
 419                 // baryon
 420                 if (kfl > 10) kfl/=10;
 421                 if (kfl > 10) kfl/=10;
 422
 423                 Int_t ipa = iparticle->GetFirstMother()-1;
 424                 Int_t kfMo = 0;
 425
 426                 if (ipa > -1) {
 427                     TParticle *  mother = (TParticle *) fParticles->At(ipa);
 428                     kfMo = TMath::Abs(mother->GetPdgCode());
 429                 }
 430                 // What to keep in Stack?
 431                 Bool_t flavorOK = kFALSE;
 432                 Bool_t selectOK = kFALSE;
 433                 if (fFeedDownOpt) {
 434                     if (kfl >= fFlavorSelect) flavorOK = kTRUE;
 435                 } else {
 436                     if (kfl > fFlavorSelect) {
 437                         nc = -1;
 438                         break;
 439                     }
 440                     if (kfl == fFlavorSelect) flavorOK = kTRUE;
 441                 }
 442                 switch (fStackFillOpt) {
 443                 case kFlavorSelection:
 444                     selectOK = kTRUE;
 445                     break;
 446                 case kParentSelection:
 447                     if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
 448                     break;
 449                 }
 450                 if (flavorOK && selectOK) {
 451 //
 452 // Heavy flavor hadron or quark
 453 //
 454 // Kinematic seletion on final state heavy flavor mesons
 455                     if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
 456                     {
 457                       continue;
 458                     }
 459                     pSelected[i] = 1;
 460                     if (ParentSelected(kf)) ++nParents; // Update parent count
 461 //                  printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
 462                 } else {
 463 // Kinematic seletion on decay products
 464                     if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
 465                         && !KinematicSelection(iparticle, 1))
 466                     {
 467                       continue;
 468                     }
 469 //
 470 // Decay products
 471 // Select if mother was selected and is not tracked
 472
 473                     if (pSelected[ipa] &&
 474                         !trackIt[ipa]  &&     // mother will be  tracked ?
 475                         kfMo !=  5 &&         // mother is b-quark, don't store fragments
 476                         kfMo !=  4 &&         // mother is c-quark, don't store fragments
 477                         kf   != 92)           // don't store string
 478                     {
 479 //
 480 // Semi-stable or de-selected: diselect decay products:
 481 //
 482 //
 483                         if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
 484                         {
 485                             Int_t ipF = iparticle->GetFirstDaughter();
 486                             Int_t ipL = iparticle->GetLastDaughter();
 487                             if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
 488                         }
 489 //                      printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
 490                         pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
 491                     }
 492                 }
 493                 if (pSelected[i] == -1) pSelected[i] = 0;
 494                 if (!pSelected[i]) continue;
 495                 // Count quarks only if you did not include fragmentation
 496                 if (fFragmentation && kf <= 10) continue;
 497                 nc++;
 498 // Decision on tracking
 499                 trackIt[i] = 0;
 500 //
 501 // Track final state particle
 502                 if (ks == 1) trackIt[i] = 1;
 503 // Track semi-stable particles
 504                 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
 505 // Track particles selected by process if undecayed.
 506                 if (fForceDecay == kNoDecay) {
 507                     if (ParentSelected(kf)) trackIt[i] = 1;
 508                 } else {
 509                     if (ParentSelected(kf)) trackIt[i] = 0;
 510                 }
 511                 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
 512 //
 513 //
 514
 515             } // particle selection loop
 516             if (nc > 0) {
 517                 for (i = 0; i<np; i++) {
 518                     if (!pSelected[i]) continue;
 519                     TParticle *  iparticle = (TParticle *) fParticles->At(i);
 520                     kf = CheckPDGCode(iparticle->GetPdgCode());
 521                     Int_t ks = iparticle->GetStatusCode();
 522                     p[0] = iparticle->Px();
 523                     p[1] = iparticle->Py();
 524                     p[2] = iparticle->Pz();
 525                     p[3] = iparticle->Energy();
 526
 527                     origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 528                     origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 529                     origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 530
 531                     Float_t tof   = kconv*iparticle->T();
 532                     Int_t ipa     = iparticle->GetFirstMother()-1;
 533                     Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
 534
 535                     PushTrack(fTrackIt*trackIt[i], iparent, kf,
 536                               p[0], p[1], p[2], p[3],
 537                               origin[0], origin[1], origin[2], tof,
 538                               polar[0], polar[1], polar[2],
 539                               kPPrimary, nt, 1., ks);
 540                     pParent[i] = nt;
 541                     KeepTrack(nt);
 542                 } //  PushTrack loop
 543             }
 544         } else {
 545             nc = GenerateMB();
 546         } // mb ?
 547
 548         if (pParent)   delete[] pParent;
 549         if (pSelected) delete[] pSelected;
 550         if (trackIt)   delete[] trackIt;
 551
 552         if (nc > 0) {
 553           switch (fCountMode) {
 554           case kCountAll:
 555             // printf(" Count all \n");
 556             jev += nc;
 557             break;
 558           case kCountParents:
 559             // printf(" Count parents \n");
 560             jev += nParents;
 561             break;
 562           case kCountTrackables:
 563             // printf(" Count trackable \n");
 564             jev += nTkbles;
 565             break;
 566           }
 567             if (jev >= fNpart || fNpart == -1) {
 568                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 569                 printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
 570
 571                 fQ  += fPythia->GetVINT(51);
 572                 fX1 += fPythia->GetVINT(41);
 573                 fX2 += fPythia->GetVINT(42);
 574                 fTrialsRun += fTrials;
 575                 fNev++;
 576                 MakeHeader();
 577                 break;
 578             }
 579         }
 580     } // event loop
 581     SetHighWaterMark(nt);
 582 //  adjust weight due to kinematic selection
 583 //    AdjustWeights();
 584 //  get cross-section
 585     fXsection=fPythia->GetPARI(1);
 586 }
 587
 588 Int_t  AliGenPythia::GenerateMB()
 589 {
 590 //
 591 // Min Bias selection and other global selections
 592 //
 593     Int_t i, kf, nt, iparent;
 594     Int_t nc = 0;
 595     Float_t p[4];
 596     Float_t polar[3]   =   {0,0,0};
 597     Float_t origin[3]  =   {0,0,0};
 598 //  converts from mm/c to s
 599     const Float_t kconv=0.001/2.999792458e8;
 600
 601
 602
 603     Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
 604
 605
 606     Int_t* pParent = new Int_t[np];
 607     for (i=0; i< np; i++) pParent[i] = -1;
 608     if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
 609         TParticle* jet1 = (TParticle *) fParticles->At(6);
 610         TParticle* jet2 = (TParticle *) fParticles->At(7);
 611         if (!CheckTrigger(jet1, jet2)) return 0;
 612     }
 613
 614     for (i = 0; i<np; i++) {
 615         Int_t trackIt = 0;
 616         TParticle *  iparticle = (TParticle *) fParticles->At(i);
 617         kf = CheckPDGCode(iparticle->GetPdgCode());
 618         Int_t ks = iparticle->GetStatusCode();
 619         Int_t km = iparticle->GetFirstMother();
 620         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
 621             (ks != 1) ||
 622             (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
 623             nc++;
 624             if (ks == 1) trackIt = 1;
 625             Int_t ipa = iparticle->GetFirstMother()-1;
 626
 627             iparent = (ipa > -1) ? pParent[ipa] : -1;
 628
 629 //
 630 // store track information
 631             p[0] = iparticle->Px();
 632             p[1] = iparticle->Py();
 633             p[2] = iparticle->Pz();
 634             p[3] = iparticle->Energy();
 635
 636
 637             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 638             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 639             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 640
 641             Float_t tof=kconv*iparticle->T();
 642
 643             PushTrack(fTrackIt*trackIt, iparent, kf,
 644                       p[0], p[1], p[2], p[3],
 645                       origin[0], origin[1], origin[2], tof,
 646                       polar[0], polar[1], polar[2],
 647                       kPPrimary, nt, 1., ks);
 648             //
 649             // Special Treatment to store color-flow
 650             //
 651             if (ks == 3 || ks == 13 || ks == 14) {
 652                 TParticle* particle = 0;
 653                 if (fStack) {
 654                     particle = fStack->Particle(nt);
 655                 } else {
 656                     particle = gAlice->Stack()->Particle(nt);
 657                 }
 658                 particle->SetFirstDaughter(fPythia->GetK(2, i));
 659                 particle->SetLastDaughter(fPythia->GetK(3, i));
 660             }
 661
 662             KeepTrack(nt);
 663             pParent[i] = nt;
 664         } // select particle
 665     } // particle loop
 666
 667     if (pParent) delete[] pParent;
 668
 669     printf("\n I've put %i particles on the stack \n",nc);
 670     return nc;
 671 }
 672
 673
 674 void AliGenPythia::FinishRun()
 675 {
 676 // Print x-section summary
 677     fPythia->Pystat(1);
 678     fQ  /= fNev;
 679     fX1 /= fNev;
 680     fX2 /= fNev;
 681     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
 682     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
 683
 684
 685 }
 686
 687 void AliGenPythia::AdjustWeights()
 688 {
 689 // Adjust the weights after generation of all events
 690 //
 691     if (gAlice) {
 692         TParticle *part;
 693         Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
 694         for (Int_t i=0; i<ntrack; i++) {
 695             part= gAlice->GetMCApp()->Particle(i);
 696             part->SetWeight(part->GetWeight()*fKineBias);
 697         }
 698     }
 699 }
 700
 701 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
 702 {
 703 // Treat protons as inside nuclei with mass numbers a1 and a2
 704
 705     fAProjectile = a1;
 706     fATarget     = a2;
 707     fSetNuclei   = kTRUE;
 708 }
 709
 710
 711 void AliGenPythia::MakeHeader()
 712 {
 713 // Builds the event header, to be called after each event
 714     if (fHeader) delete fHeader;
 715     fHeader = new AliGenPythiaEventHeader("Pythia");
 716 //
 717 // Event type
 718     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
 719 //
 720 // Number of trials
 721     ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
 722 //
 723 // Event Vertex
 724     fHeader->SetPrimaryVertex(fVertex);
 725 //
 726 // Jets that have triggered
 727     if (fProcess == kPyJets)
 728     {
 729         Int_t ntrig, njet;
 730         Float_t jets[4][10];
 731         GetJets(njet, ntrig, jets);
 732
 733         for (Int_t i = 0; i < ntrig; i++) {
 734             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
 735                                                         jets[3][i]);
 736         }
 737     }
 738 //
 739 // Copy relevant information from external header, if present.
 740 //
 741     Float_t uqJet[4];
 742
 743     if (fRL) {
 744         AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
 745         for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
 746         {
 747             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
 748
 749
 750             exHeader->TriggerJet(i, uqJet);
 751             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
 752         }
 753     }
 754 //
 755 // Store quenching parameters
 756 //
 757     if (fQuench){
 758         Double_t z[4];
 759         Double_t xp, yp;
 760
 761         fPythia->GetQuenchingParameters(xp, yp, z);
 762
 763         ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
 764         ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
 765     }
 766
 767 //
 768 //  Pass header to RunLoader
 769 //
 770     AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(fHeader);
 771 }
 772
 773
 774 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
 775 {
 776 // Check the kinematic trigger condition
 777 //
 778     Double_t eta[2];
 779     eta[0] = jet1->Eta();
 780     eta[1] = jet2->Eta();
 781     Double_t phi[2];
 782     phi[0] = jet1->Phi();
 783     phi[1] = jet2->Phi();
 784     Int_t    pdg[2];
 785     pdg[0] = jet1->GetPdgCode();
 786     pdg[1] = jet2->GetPdgCode();
 787     Bool_t   triggered = kFALSE;
 788
 789     if (fProcess == kPyJets) {
 790         Int_t njets = 0;
 791         Int_t ntrig = 0;
 792         Float_t jets[4][10];
 793 //
 794 // Use Pythia clustering on parton level to determine jet axis
 795 //
 796         GetJets(njets, ntrig, jets);
 797
 798         if (ntrig) triggered = kTRUE;
 799 //
 800     } else {
 801         Int_t ij = 0;
 802         Int_t ig = 1;
 803         if (pdg[0] == kGamma) {
 804             ij = 1;
 805             ig = 0;
 806         }
 807         //Check eta range first...
 808         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
 809             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
 810         {
 811             //Eta is okay, now check phi range
 812             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
 813                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
 814             {
 815                 triggered = kTRUE;
 816             }
 817         }
 818     }
 819     return triggered;
 820 }
 821
 822 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
 823 {
 824 // Assignment operator
 825     rhs.Copy(*this);
 826     return *this;
 827 }
 828
 829 void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
 830 {
 831 //
 832 // Load event into Pythia Common Block
 833 //
 834
 835     Int_t npart = stack -> GetNprimary();
 836     Int_t n0 = 0;
 837
 838     if (!flag) {
 839         (fPythia->GetPyjets())->N = npart;
 840     } else {
 841         n0 = (fPythia->GetPyjets())->N;
 842         (fPythia->GetPyjets())->N = n0 + npart;
 843     }
 844
 845
 846     for (Int_t part = 0; part < npart; part++) {
 847         TParticle *MPart = stack->Particle(part);
 848
 849         Int_t kf     =  MPart->GetPdgCode();
 850         Int_t ks     =  MPart->GetStatusCode();
 851         Int_t idf    =  MPart->GetFirstDaughter();
 852         Int_t idl    =  MPart->GetLastDaughter();
 853
 854         if (reHadr) {
 855             if (ks == 11 || ks == 12) {
 856                 ks  -= 10;
 857                 idf  = -1;
 858                 idl  = -1;
 859             }
 860         }
 861
 862         Float_t px = MPart->Px();
 863         Float_t py = MPart->Py();
 864         Float_t pz = MPart->Pz();
 865         Float_t e  = MPart->Energy();
 866         Float_t m  = MPart->GetCalcMass();
 867
 868
 869         (fPythia->GetPyjets())->P[0][part+n0] = px;
 870         (fPythia->GetPyjets())->P[1][part+n0] = py;
 871         (fPythia->GetPyjets())->P[2][part+n0] = pz;
 872         (fPythia->GetPyjets())->P[3][part+n0] = e;
 873         (fPythia->GetPyjets())->P[4][part+n0] = m;
 874
 875         (fPythia->GetPyjets())->K[1][part+n0] = kf;
 876         (fPythia->GetPyjets())->K[0][part+n0] = ks;
 877         (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
 878         (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
 879         (fPythia->GetPyjets())->K[2][part+n0] = MPart->GetFirstMother() + 1;
 880     }
 881 }
 882
 883
 884 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
 885 {
 886 //
 887 //  Calls the Pythia jet finding algorithm to find jets in the current event
 888 //
 889 //
 890 //
 891 //  Save jets
 892     Int_t n     = fPythia->GetN();
 893
 894 //
 895 //  Run Jet Finder
 896     fPythia->Pycell(njets);
 897     Int_t i;
 898     for (i = 0; i < njets; i++) {
 899         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
 900         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
 901         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
 902         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
 903
 904         jets[0][i] = px;
 905         jets[1][i] = py;
 906         jets[2][i] = pz;
 907         jets[3][i] = e;
 908     }
 909 }
 910
 911
 912
 913 void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
 914 {
 915 //
 916 //  Calls the Pythia clustering algorithm to find jets in the current event
 917 //
 918     Int_t n     = fPythia->GetN();
 919     nJets       = 0;
 920     nJetsTrig   = 0;
 921     if (fJetReconstruction == kCluster) {
 922 //
 923 //  Configure cluster algorithm
 924 //
 925         fPythia->SetPARU(43, 2.);
 926         fPythia->SetMSTU(41, 1);
 927 //
 928 //  Call cluster algorithm
 929 //
 930         fPythia->Pyclus(nJets);
 931 //
 932 //  Loading jets from common block
 933 //
 934     } else {
 935
 936 //
 937 //  Run Jet Finder
 938         fPythia->Pycell(nJets);
 939     }
 940
 941     Int_t i;
 942     for (i = 0; i < nJets; i++) {
 943         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
 944         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
 945         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
 946         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
 947         Float_t pt    = TMath::Sqrt(px * px + py * py);
 948         Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);
 949         Float_t theta = TMath::ATan2(pt,pz);
 950         Float_t et    = e * TMath::Sin(theta);
 951         Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
 952
 953         if (
 954             eta > fEtaMinJet && eta < fEtaMaxJet &&
 955             phi > fPhiMinJet && phi < fPhiMaxJet &&
 956             et  > fEtMinJet  && et  < fEtMaxJet
 957             )
 958         {
 959             jets[0][nJetsTrig] = px;
 960             jets[1][nJetsTrig] = py;
 961             jets[2][nJetsTrig] = pz;
 962             jets[3][nJetsTrig] = e;
 963             nJetsTrig++;
 964 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
 965         } else {
 966 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
 967         }
 968     }
 969 }
 970
 971
 972 #ifdef never
 973 void AliGenPythia::Streamer(TBuffer &R__b)
 974 {
 975    // Stream an object of class AliGenPythia.
 976
 977    if (R__b.IsReading()) {
 978       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
 979       AliGenerator::Streamer(R__b);
 980       R__b >> (Int_t&)fProcess;
 981       R__b >> (Int_t&)fStrucFunc;
 982       R__b >> (Int_t&)fForceDecay;
 983       R__b >> fEnergyCMS;
 984       R__b >> fKineBias;
 985       R__b >> fTrials;
 986       fParentSelect.Streamer(R__b);
 987       fChildSelect.Streamer(R__b);
 988       R__b >> fXsection;
 989 //      (AliPythia::Instance())->Streamer(R__b);
 990       R__b >> fPtHardMin;
 991       R__b >> fPtHardMax;
 992 //      if (fDecayer) fDecayer->Streamer(R__b);
 993    } else {
 994       R__b.WriteVersion(AliGenPythia::IsA());
 995       AliGenerator::Streamer(R__b);
 996       R__b << (Int_t)fProcess;
 997       R__b << (Int_t)fStrucFunc;
 998       R__b << (Int_t)fForceDecay;
 999       R__b << fEnergyCMS;
1000       R__b << fKineBias;
1001       R__b << fTrials;
1002       fParentSelect.Streamer(R__b);
1003       fChildSelect.Streamer(R__b);
1004       R__b << fXsection;
1005 //      R__b << fPythia;
1006       R__b << fPtHardMin;
1007       R__b << fPtHardMax;
1008       //     fDecayer->Streamer(R__b);
1009    }
1010 }
1011 #endif
1012
1013