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 <TClonesArray.h>
  29 #include <TDatabasePDG.h>
  30 #include <TParticle.h>
  31 #include <TPDGCode.h>
  32 #include <TObjArray.h>
  33 #include <TSystem.h>
  34 #include <TTree.h>
  35 #include "AliConst.h"
  36 #include "AliDecayerPythia.h"
  37 #include "AliGenPythia.h"
  38 #include "AliFastGlauber.h"
  39 #include "AliHeader.h"
  40 #include "AliGenPythiaEventHeader.h"
  41 #include "AliPythia.h"
  42 #include "AliPythiaRndm.h"
  43 #include "AliRun.h"
  44 #include "AliStack.h"
  45 #include "AliRunLoader.h"
  46 #include "AliMC.h"
  47 #include "AliLog.h"
  48 #include "PyquenCommon.h"
  49 #include "AliLog.h"
  50
  51 ClassImp(AliGenPythia)
  52
  53
  54 AliGenPythia::AliGenPythia():
  55     AliGenMC(),
  56     fProcess(kPyCharm),
  57     fItune(-1),
  58     fStrucFunc(kCTEQ5L),
  59     fKineBias(0.),
  60     fTrials(0),
  61     fTrialsRun(0),
  62     fQ(0.),
  63     fX1(0.),
  64     fX2(0.),
  65     fEventTime(0.),
  66     fInteractionRate(0.),
  67     fTimeWindow(0.),
  68     fCurSubEvent(0),
  69     fEventsTime(0),
  70     fNev(0),
  71     fFlavorSelect(0),
  72     fXsection(0.),
  73     fPythia(0),
  74     fPtHardMin(0.),
  75     fPtHardMax(1.e4),
  76     fYHardMin(-1.e10),
  77     fYHardMax(1.e10),
  78     fGinit(1),
  79     fGfinal(1),
  80     fHadronisation(1),
  81     fPatchOmegaDalitz(0),
  82     fNpartons(0),
  83     fReadFromFile(0),
  84     fQuench(0),
  85     fQhat(0.),
  86     fLength(0.),
  87     fImpact(0.),
  88     fPtKick(1.),
  89     fFullEvent(kTRUE),
  90     fDecayer(new AliDecayerPythia()),
  91     fDebugEventFirst(-1),
  92     fDebugEventLast(-1),
  93     fEtMinJet(0.),
  94     fEtMaxJet(1.e4),
  95     fEtaMinJet(-20.),
  96     fEtaMaxJet(20.),
  97     fPhiMinJet(0.),
  98     fPhiMaxJet(2.* TMath::Pi()),
  99     fJetReconstruction(kCell),
 100     fEtaMinGamma(-20.),
 101     fEtaMaxGamma(20.),
 102     fPhiMinGamma(0.),
 103     fPhiMaxGamma(2. * TMath::Pi()),
 104     fPycellEtaMax(2.),
 105     fPycellNEta(274),
 106     fPycellNPhi(432),
 107     fPycellThreshold(0.),
 108     fPycellEtSeed(4.),
 109     fPycellMinEtJet(10.),
 110     fPycellMaxRadius(1.),
 111     fStackFillOpt(kFlavorSelection),
 112     fFeedDownOpt(kTRUE),
 113     fFragmentation(kTRUE),
 114     fSetNuclei(kFALSE),
 115     fNewMIS(kFALSE),
 116     fHFoff(kFALSE),
 117     fNucPdf(0),
 118     fTriggerParticle(0),
 119     fTriggerEta(0.9),
 120     fTriggerMultiplicity(0),
 121     fTriggerMultiplicityEta(0),
 122     fTriggerMultiplicityPtMin(0),
 123     fCountMode(kCountAll),
 124     fHeader(0),
 125     fRL(0),
 126     fkFileName(0),
 127     fFragPhotonInCalo(kFALSE),
 128     fPi0InCalo(kFALSE) ,
 129     fPhotonInCalo(kFALSE),
 130     fEleInEMCAL(kFALSE),
 131     fCheckEMCAL(kFALSE),
 132     fCheckPHOS(kFALSE),
 133     fCheckPHOSeta(kFALSE),
 134     fFragPhotonOrPi0MinPt(0),
 135     fPhotonMinPt(0),
 136     fElectronMinPt(0),
 137     fPHOSMinPhi(219.),
 138     fPHOSMaxPhi(321.),
 139     fPHOSEta(0.13),
 140     fEMCALMinPhi(79.),
 141     fEMCALMaxPhi(191.),
 142     fEMCALEta(0.71),
 143     fkTuneForDiff(0),
 144     fProcDiff(0)
 145 {
 146 // Default Constructor
 147   fEnergyCMS = 5500.;
 148   if (!AliPythiaRndm::GetPythiaRandom())
 149       AliPythiaRndm::SetPythiaRandom(GetRandom());
 150 }
 151
 152 AliGenPythia::AliGenPythia(Int_t npart)
 153     :AliGenMC(npart),
 154      fProcess(kPyCharm),
 155      fItune(-1),
 156      fStrucFunc(kCTEQ5L),
 157      fKineBias(0.),
 158      fTrials(0),
 159      fTrialsRun(0),
 160      fQ(0.),
 161      fX1(0.),
 162      fX2(0.),
 163      fEventTime(0.),
 164      fInteractionRate(0.),
 165      fTimeWindow(0.),
 166      fCurSubEvent(0),
 167      fEventsTime(0),
 168      fNev(0),
 169      fFlavorSelect(0),
 170      fXsection(0.),
 171      fPythia(0),
 172      fPtHardMin(0.),
 173      fPtHardMax(1.e4),
 174      fYHardMin(-1.e10),
 175      fYHardMax(1.e10),
 176      fGinit(kTRUE),
 177      fGfinal(kTRUE),
 178      fHadronisation(kTRUE),
 179      fPatchOmegaDalitz(0),
 180      fNpartons(0),
 181      fReadFromFile(kFALSE),
 182      fQuench(kFALSE),
 183      fQhat(0.),
 184      fLength(0.),
 185      fImpact(0.),
 186      fPtKick(1.),
 187      fFullEvent(kTRUE),
 188      fDecayer(new AliDecayerPythia()),
 189      fDebugEventFirst(-1),
 190      fDebugEventLast(-1),
 191      fEtMinJet(0.),
 192      fEtMaxJet(1.e4),
 193      fEtaMinJet(-20.),
 194      fEtaMaxJet(20.),
 195      fPhiMinJet(0.),
 196      fPhiMaxJet(2.* TMath::Pi()),
 197      fJetReconstruction(kCell),
 198      fEtaMinGamma(-20.),
 199      fEtaMaxGamma(20.),
 200      fPhiMinGamma(0.),
 201      fPhiMaxGamma(2. * TMath::Pi()),
 202      fPycellEtaMax(2.),
 203      fPycellNEta(274),
 204      fPycellNPhi(432),
 205      fPycellThreshold(0.),
 206      fPycellEtSeed(4.),
 207      fPycellMinEtJet(10.),
 208      fPycellMaxRadius(1.),
 209      fStackFillOpt(kFlavorSelection),
 210      fFeedDownOpt(kTRUE),
 211      fFragmentation(kTRUE),
 212      fSetNuclei(kFALSE),
 213      fNewMIS(kFALSE),
 214      fHFoff(kFALSE),
 215      fNucPdf(0),
 216      fTriggerParticle(0),
 217      fTriggerEta(0.9),
 218      fTriggerMultiplicity(0),
 219      fTriggerMultiplicityEta(0),
 220      fTriggerMultiplicityPtMin(0),
 221      fCountMode(kCountAll),
 222      fHeader(0),
 223      fRL(0),
 224      fkFileName(0),
 225      fFragPhotonInCalo(kFALSE),
 226      fPi0InCalo(kFALSE) ,
 227      fPhotonInCalo(kFALSE),
 228      fEleInEMCAL(kFALSE),
 229      fCheckEMCAL(kFALSE),
 230      fCheckPHOS(kFALSE),
 231      fCheckPHOSeta(kFALSE),
 232      fFragPhotonOrPi0MinPt(0),
 233      fPhotonMinPt(0),
 234      fElectronMinPt(0),
 235      fPHOSMinPhi(219.),
 236      fPHOSMaxPhi(321.),
 237      fPHOSEta(0.13),
 238      fEMCALMinPhi(79.),
 239      fEMCALMaxPhi(191.),
 240      fEMCALEta(0.71),
 241      fkTuneForDiff(0),
 242      fProcDiff(0)
 243 {
 244 // default charm production at 5. 5 TeV
 245 // semimuonic decay
 246 // structure function GRVHO
 247 //
 248     fEnergyCMS = 5500.;
 249     fName = "Pythia";
 250     fTitle= "Particle Generator using PYTHIA";
 251     SetForceDecay();
 252     // Set random number generator
 253     if (!AliPythiaRndm::GetPythiaRandom())
 254       AliPythiaRndm::SetPythiaRandom(GetRandom());
 255  }
 256
 257 AliGenPythia::~AliGenPythia()
 258 {
 259 // Destructor
 260   if(fEventsTime) delete fEventsTime;
 261 }
 262
 263 void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
 264 {
 265 // Generate pileup using user specified rate
 266     fInteractionRate = rate;
 267     fTimeWindow = timewindow;
 268     GeneratePileup();
 269 }
 270
 271 void AliGenPythia::GeneratePileup()
 272 {
 273 // Generate sub events time for pileup
 274     fEventsTime = 0;
 275     if(fInteractionRate == 0.) {
 276       Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
 277       return;
 278     }
 279
 280     Int_t npart = NumberParticles();
 281     if(npart < 0) {
 282       Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
 283       return;
 284     }
 285
 286     if(fEventsTime) delete fEventsTime;
 287     fEventsTime = new TArrayF(npart);
 288     TArrayF &array = *fEventsTime;
 289     for(Int_t ipart = 0; ipart < npart; ipart++)
 290       array[ipart] = 0.;
 291
 292     Float_t eventtime = 0.;
 293     while(1)
 294       {
 295         eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
 296         if(eventtime > fTimeWindow) break;
 297         array.Set(array.GetSize()+1);
 298         array[array.GetSize()-1] = eventtime;
 299       }
 300
 301     eventtime = 0.;
 302     while(1)
 303       {
 304         eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
 305         if(TMath::Abs(eventtime) > fTimeWindow) break;
 306         array.Set(array.GetSize()+1);
 307         array[array.GetSize()-1] = eventtime;
 308       }
 309
 310     SetNumberParticles(fEventsTime->GetSize());
 311 }
 312
 313 void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
 314                                        Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
 315 {
 316 // Set pycell parameters
 317     fPycellEtaMax    =  etamax;
 318     fPycellNEta      =  neta;
 319     fPycellNPhi      =  nphi;
 320     fPycellThreshold =  thresh;
 321     fPycellEtSeed    =  etseed;
 322     fPycellMinEtJet  =  minet;
 323     fPycellMaxRadius =  r;
 324 }
 325
 326
 327
 328 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
 329 {
 330   // Set a range of event numbers, for which a table
 331   // of generated particle will be printed
 332   fDebugEventFirst = eventFirst;
 333   fDebugEventLast  = eventLast;
 334   if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
 335 }
 336
 337 void AliGenPythia::Init()
 338 {
 339 // Initialisation
 340
 341     SetMC(AliPythia::Instance());
 342     fPythia=(AliPythia*) fMCEvGen;
 343
 344 //
 345     fParentWeight=1./Float_t(fNpart);
 346 //
 347
 348
 349     fPythia->SetCKIN(3,fPtHardMin);
 350     fPythia->SetCKIN(4,fPtHardMax);
 351     fPythia->SetCKIN(7,fYHardMin);
 352     fPythia->SetCKIN(8,fYHardMax);
 353
 354     if (fAProjectile > 0 && fATarget > 0) fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);
 355     // Fragmentation?
 356     if (fFragmentation) {
 357       fPythia->SetMSTP(111,1);
 358     } else {
 359       fPythia->SetMSTP(111,0);
 360     }
 361
 362
 363 //  initial state radiation
 364     fPythia->SetMSTP(61,fGinit);
 365 //  final state radiation
 366     fPythia->SetMSTP(71,fGfinal);
 367 //  pt - kick
 368     if (fPtKick > 0.) {
 369         fPythia->SetMSTP(91,1);
 370         fPythia->SetPARP(91,fPtKick);
 371         fPythia->SetPARP(93, 4. * fPtKick);
 372     } else {
 373         fPythia->SetMSTP(91,0);
 374     }
 375
 376
 377     if (fReadFromFile) {
 378         fRL  =  AliRunLoader::Open(fkFileName, "Partons");
 379         fRL->LoadKinematics();
 380         fRL->LoadHeader();
 381     } else {
 382         fRL = 0x0;
 383     }
 384  //
 385     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc, fItune);
 386     //  Forward Paramters to the AliPythia object
 387     fDecayer->SetForceDecay(fForceDecay);
 388 // Switch off Heavy Flavors on request
 389     if (fHFoff) {
 390         // Maximum number of quark flavours used in pdf
 391         fPythia->SetMSTP(58, 3);
 392         // Maximum number of flavors that can be used in showers
 393         fPythia->SetMSTJ(45, 3);
 394         // Switch off g->QQbar splitting in decay table
 395         ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
 396     }
 397
 398     fDecayer->Init();
 399
 400
 401 //  Parent and Children Selection
 402     switch (fProcess)
 403     {
 404     case kPyOldUEQ2ordered:
 405     case kPyOldUEQ2ordered2:
 406     case kPyOldPopcorn:
 407       break;
 408     case kPyCharm:
 409     case kPyCharmUnforced:
 410     case kPyCharmPbPbMNR:
 411     case kPyCharmpPbMNR:
 412     case kPyCharmppMNR:
 413     case kPyCharmppMNRwmi:
 414         fParentSelect[0] =   411;
 415         fParentSelect[1] =   421;
 416         fParentSelect[2] =   431;
 417         fParentSelect[3] =  4122;
 418         fParentSelect[4] =  4232;
 419         fParentSelect[5] =  4132;
 420         fParentSelect[6] =  4332;
 421         fFlavorSelect    =  4;
 422         break;
 423     case kPyD0PbPbMNR:
 424     case kPyD0pPbMNR:
 425     case kPyD0ppMNR:
 426         fParentSelect[0] =   421;
 427         fFlavorSelect    =   4;
 428         break;
 429     case kPyDPlusPbPbMNR:
 430     case kPyDPluspPbMNR:
 431     case kPyDPlusppMNR:
 432         fParentSelect[0] =   411;
 433         fFlavorSelect    =   4;
 434         break;
 435     case kPyDPlusStrangePbPbMNR:
 436     case kPyDPlusStrangepPbMNR:
 437     case kPyDPlusStrangeppMNR:
 438         fParentSelect[0] =   431;
 439         fFlavorSelect    =   4;
 440         break;
 441     case kPyLambdacppMNR:
 442         fParentSelect[0] =  4122;
 443         fFlavorSelect    =   4;
 444         break;
 445     case kPyBeauty:
 446     case kPyBeautyJets:
 447     case kPyBeautyPbPbMNR:
 448     case kPyBeautypPbMNR:
 449     case kPyBeautyppMNR:
 450     case kPyBeautyppMNRwmi:
 451         fParentSelect[0]=  511;
 452         fParentSelect[1]=  521;
 453         fParentSelect[2]=  531;
 454         fParentSelect[3]= 5122;
 455         fParentSelect[4]= 5132;
 456         fParentSelect[5]= 5232;
 457         fParentSelect[6]= 5332;
 458         fFlavorSelect   = 5;
 459         break;
 460     case kPyBeautyUnforced:
 461         fParentSelect[0] =  511;
 462         fParentSelect[1] =  521;
 463         fParentSelect[2] =  531;
 464         fParentSelect[3] = 5122;
 465         fParentSelect[4] = 5132;
 466         fParentSelect[5] = 5232;
 467         fParentSelect[6] = 5332;
 468         fFlavorSelect    = 5;
 469         break;
 470     case kPyJpsiChi:
 471     case kPyJpsi:
 472         fParentSelect[0] = 443;
 473         break;
 474     case kPyMbDefault:
 475     case kPyMbAtlasTuneMC09:
 476     case kPyMb:
 477     case kPyMbWithDirectPhoton:
 478     case kPyMbNonDiffr:
 479     case kPyMbMSEL1:
 480     case kPyJets:
 481     case kPyDirectGamma:
 482     case kPyLhwgMb:
 483         break;
 484     case kPyW:
 485     case kPyZ:
 486         break;
 487     }
 488 //
 489 //
 490 //  JetFinder for Trigger
 491 //
 492 //  Configure detector (EMCAL like)
 493 //
 494     fPythia->SetPARU(51, fPycellEtaMax);
 495     fPythia->SetMSTU(51, fPycellNEta);
 496     fPythia->SetMSTU(52, fPycellNPhi);
 497 //
 498 //  Configure Jet Finder
 499 //
 500     fPythia->SetPARU(58,  fPycellThreshold);
 501     fPythia->SetPARU(52,  fPycellEtSeed);
 502     fPythia->SetPARU(53,  fPycellMinEtJet);
 503     fPythia->SetPARU(54,  fPycellMaxRadius);
 504     fPythia->SetMSTU(54,  2);
 505 //
 506 //  This counts the total number of calls to Pyevnt() per run.
 507     fTrialsRun = 0;
 508     fQ         = 0.;
 509     fX1        = 0.;
 510     fX2        = 0.;
 511     fNev       = 0 ;
 512 //
 513 //
 514 //
 515     AliGenMC::Init();
 516 //
 517 //
 518 //
 519     if (fSetNuclei) {
 520         fDyBoost = 0;
 521         Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
 522     }
 523
 524     fPythia->SetPARJ(200, 0.0);
 525     fPythia->SetPARJ(199, 0.0);
 526     fPythia->SetPARJ(198, 0.0);
 527     fPythia->SetPARJ(197, 0.0);
 528
 529     if (fQuench == 1) {
 530         fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
 531     }
 532
 533     if (fQuench == 3) {
 534         // Nestor's change of the splittings
 535         fPythia->SetPARJ(200, 0.8);
 536         fPythia->SetMSTJ(41, 1);  // QCD radiation only
 537         fPythia->SetMSTJ(42, 2);  // angular ordering
 538         fPythia->SetMSTJ(44, 2);  // option to run alpha_s
 539         fPythia->SetMSTJ(47, 0);  // No correction back to hard scattering element
 540         fPythia->SetMSTJ(50, 0);  // No coherence in first branching
 541         fPythia->SetPARJ(82, 1.); // Cut off for parton showers
 542     } else if (fQuench == 4) {
 543         // Armesto-Cunqueiro-Salgado change of the splittings.
 544         AliFastGlauber* glauber = AliFastGlauber::Instance();
 545         glauber->Init(2);
 546         //read and store transverse almonds corresponding to differnt
 547         //impact parameters.
 548         glauber->SetCentralityClass(0.,0.1);
 549         fPythia->SetPARJ(200, 1.);
 550         fPythia->SetPARJ(198, fQhat);
 551         fPythia->SetPARJ(199, fLength);
 552         fPythia->SetMSTJ(42, 2);  // angular ordering
 553         fPythia->SetMSTJ(44, 2);  // option to run alpha_s
 554         fPythia->SetPARJ(82, 1.); // Cut off for parton showers
 555     }
 556 }
 557
 558 void AliGenPythia::Generate()
 559 {
 560 // Generate one event
 561     if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
 562     fDecayer->ForceDecay();
 563
 564     Double_t polar[3]   =   {0,0,0};
 565     Double_t origin[3]  =   {0,0,0};
 566     Double_t p[4];
 567 //  converts from mm/c to s
 568     const Float_t kconv=0.001/2.999792458e8;
 569 //
 570     Int_t nt=0;
 571     Int_t jev=0;
 572     Int_t j, kf;
 573     fTrials=0;
 574     fEventTime = 0.;
 575
 576
 577
 578     //  Set collision vertex position
 579     if (fVertexSmear == kPerEvent) Vertex();
 580
 581 //  event loop
 582     while(1)
 583     {
 584 //
 585 // Produce event
 586 //
 587 //
 588 // Switch hadronisation off
 589 //
 590         fPythia->SetMSTJ(1, 0);
 591
 592         if (fQuench ==4){
 593             Double_t bimp;
 594             // Quenching comes through medium-modified splitting functions.
 595             AliFastGlauber::Instance()->GetRandomBHard(bimp);
 596             fPythia->SetPARJ(197, bimp);
 597             fImpact = bimp;
 598             fPythia->Qpygin0();
 599         }
 600 //
 601 // Either produce new event or read partons from file
 602 //
 603         if (!fReadFromFile) {
 604             if (!fNewMIS) {
 605                 fPythia->Pyevnt();
 606             } else {
 607                 fPythia->Pyevnw();
 608             }
 609             fNpartons = fPythia->GetN();
 610         } else {
 611             printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
 612             fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
 613             fPythia->SetN(0);
 614             LoadEvent(fRL->Stack(), 0 , 1);
 615             fPythia->Pyedit(21);
 616         }
 617
 618 //
 619 //  Run quenching routine
 620 //
 621         if (fQuench == 1) {
 622             fPythia->Quench();
 623         } else if (fQuench == 2){
 624             fPythia->Pyquen(208., 0, 0.);
 625         } else if (fQuench == 3) {
 626             // Quenching is via multiplicative correction of the splittings
 627         }
 628
 629 //
 630 // Switch hadronisation on
 631 //
 632         if (fHadronisation) {
 633             fPythia->SetMSTJ(1, 1);
 634 //
 635 // .. and perform hadronisation
 636 //      printf("Calling hadronisation %d\n", fPythia->GetN());
 637
 638             if (fPatchOmegaDalitz) {
 639               fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
 640               fPythia->Pyexec();
 641               fPythia->DalitzDecays();
 642               fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
 643             }
 644             fPythia->Pyexec();
 645         }
 646         fTrials++;
 647         fPythia->ImportParticles(&fParticles,"All");
 648
 649         if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
 650         if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
 651
 652 //
 653 //
 654 //
 655         Int_t i;
 656
 657         fNprimaries = 0;
 658         Int_t np = fParticles.GetEntriesFast();
 659
 660         if (np == 0) continue;
 661 //
 662
 663 //
 664         Int_t* pParent   = new Int_t[np];
 665         Int_t* pSelected = new Int_t[np];
 666         Int_t* trackIt   = new Int_t[np];
 667         for (i = 0; i < np; i++) {
 668             pParent[i]   = -1;
 669             pSelected[i] =  0;
 670             trackIt[i]   =  0;
 671         }
 672
 673         Int_t nc = 0;        // Total n. of selected particles
 674         Int_t nParents = 0;  // Selected parents
 675         Int_t nTkbles = 0;   // Trackable particles
 676         if (fProcess != kPyMbDefault &&
 677             fProcess != kPyMb &&
 678             fProcess != kPyMbAtlasTuneMC09 &&
 679             fProcess != kPyMbWithDirectPhoton &&
 680             fProcess != kPyJets &&
 681             fProcess != kPyDirectGamma &&
 682             fProcess != kPyMbNonDiffr  &&
 683             fProcess != kPyMbMSEL1     &&
 684             fProcess != kPyW &&
 685             fProcess != kPyZ &&
 686             fProcess != kPyCharmppMNRwmi &&
 687             fProcess != kPyBeautyppMNRwmi &&
 688             fProcess != kPyBeautyJets) {
 689
 690             for (i = 0; i < np; i++) {
 691                 TParticle* iparticle = (TParticle *) fParticles.At(i);
 692                 Int_t ks = iparticle->GetStatusCode();
 693                 kf = CheckPDGCode(iparticle->GetPdgCode());
 694 // No initial state partons
 695                 if (ks==21) continue;
 696 //
 697 // Heavy Flavor Selection
 698 //
 699                 // quark ?
 700                 kf = TMath::Abs(kf);
 701                 Int_t kfl = kf;
 702                 // Resonance
 703
 704                 if (kfl > 100000) kfl %= 100000;
 705                 if (kfl > 10000)  kfl %= 10000;
 706                 // meson ?
 707                 if  (kfl > 10) kfl/=100;
 708                 // baryon
 709                 if (kfl > 10) kfl/=10;
 710                 Int_t ipa = iparticle->GetFirstMother()-1;
 711                 Int_t kfMo = 0;
 712 //
 713 // Establish mother daughter relation between heavy quarks and mesons
 714 //
 715                 if (kf >= fFlavorSelect && kf <= 6) {
 716                     Int_t idau = iparticle->GetFirstDaughter() - 1;
 717                     if (idau > -1) {
 718                         TParticle* daughter = (TParticle *) fParticles.At(idau);
 719                         Int_t pdgD = daughter->GetPdgCode();
 720                         if (pdgD == 91 || pdgD == 92) {
 721                             Int_t jmin = daughter->GetFirstDaughter() - 1;
 722                             Int_t jmax = daughter->GetLastDaughter()  - 1;
 723                             for (Int_t jp = jmin; jp <= jmax; jp++)
 724                                 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
 725                         } // is string or cluster
 726                     } // has daughter
 727                 } // heavy quark
 728
 729
 730                 if (ipa > -1) {
 731                     TParticle *  mother = (TParticle *) fParticles.At(ipa);
 732                     kfMo = TMath::Abs(mother->GetPdgCode());
 733                 }
 734
 735                 // What to keep in Stack?
 736                 Bool_t flavorOK = kFALSE;
 737                 Bool_t selectOK = kFALSE;
 738                 if (fFeedDownOpt) {
 739                     if (kfl >= fFlavorSelect) flavorOK = kTRUE;
 740                 } else {
 741                     if (kfl > fFlavorSelect) {
 742                         nc = -1;
 743                         break;
 744                     }
 745                     if (kfl == fFlavorSelect) flavorOK = kTRUE;
 746                 }
 747                 switch (fStackFillOpt) {
 748                 case kFlavorSelection:
 749                     selectOK = kTRUE;
 750                     break;
 751                 case kParentSelection:
 752                     if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
 753                     break;
 754                 }
 755                 if (flavorOK && selectOK) {
 756 //
 757 // Heavy flavor hadron or quark
 758 //
 759 // Kinematic seletion on final state heavy flavor mesons
 760                     if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
 761                     {
 762                         continue;
 763                     }
 764                     pSelected[i] = 1;
 765                     if (ParentSelected(kf)) ++nParents; // Update parent count
 766 //                  printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
 767                 } else {
 768 // Kinematic seletion on decay products
 769                     if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
 770                         && !KinematicSelection(iparticle, 1))
 771                     {
 772                         continue;
 773                     }
 774 //
 775 // Decay products
 776 // Select if mother was selected and is not tracked
 777
 778                     if (pSelected[ipa] &&
 779                         !trackIt[ipa]  &&     // mother will be  tracked ?
 780                         kfMo !=  5 &&         // mother is b-quark, don't store fragments
 781                         kfMo !=  4 &&         // mother is c-quark, don't store fragments
 782                         kf   != 92)           // don't store string
 783                     {
 784 //
 785 // Semi-stable or de-selected: diselect decay products:
 786 //
 787 //
 788                         if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
 789                         {
 790                             Int_t ipF = iparticle->GetFirstDaughter();
 791                             Int_t ipL = iparticle->GetLastDaughter();
 792                             if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
 793                         }
 794 //                      printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
 795                         pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
 796                     }
 797                 }
 798                 if (pSelected[i] == -1) pSelected[i] = 0;
 799                 if (!pSelected[i]) continue;
 800                 // Count quarks only if you did not include fragmentation
 801                 if (fFragmentation && kf <= 10) continue;
 802
 803                 nc++;
 804 // Decision on tracking
 805                 trackIt[i] = 0;
 806 //
 807 // Track final state particle
 808                 if (ks == 1) trackIt[i] = 1;
 809 // Track semi-stable particles
 810                 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
 811 // Track particles selected by process if undecayed.
 812                 if (fForceDecay == kNoDecay) {
 813                     if (ParentSelected(kf)) trackIt[i] = 1;
 814                 } else {
 815                     if (ParentSelected(kf)) trackIt[i] = 0;
 816                 }
 817                 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
 818 //
 819 //
 820
 821             } // particle selection loop
 822             if (nc > 0) {
 823                 for (i = 0; i<np; i++) {
 824                     if (!pSelected[i]) continue;
 825                     TParticle *  iparticle = (TParticle *) fParticles.At(i);
 826                     kf = CheckPDGCode(iparticle->GetPdgCode());
 827                     Int_t ks = iparticle->GetStatusCode();
 828                     p[0] = iparticle->Px();
 829                     p[1] = iparticle->Py();
 830                     p[2] = iparticle->Pz();
 831                     p[3] = iparticle->Energy();
 832
 833                     origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 834                     origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 835                     origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 836
 837                     Float_t tof   = kconv*iparticle->T();
 838                     Int_t ipa     = iparticle->GetFirstMother()-1;
 839                     Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
 840
 841                     PushTrack(fTrackIt*trackIt[i], iparent, kf,
 842                               p[0], p[1], p[2], p[3],
 843                               origin[0], origin[1], origin[2], tof,
 844                               polar[0], polar[1], polar[2],
 845                               kPPrimary, nt, 1., ks);
 846                     pParent[i] = nt;
 847                     KeepTrack(nt);
 848                     fNprimaries++;
 849                 } //  PushTrack loop
 850             }
 851         } else {
 852             nc = GenerateMB();
 853         } // mb ?
 854
 855         GetSubEventTime();
 856
 857         delete[] pParent;
 858         delete[] pSelected;
 859         delete[] trackIt;
 860
 861         if (nc > 0) {
 862           switch (fCountMode) {
 863           case kCountAll:
 864             // printf(" Count all \n");
 865             jev += nc;
 866             break;
 867           case kCountParents:
 868             // printf(" Count parents \n");
 869             jev += nParents;
 870             break;
 871           case kCountTrackables:
 872             // printf(" Count trackable \n");
 873             jev += nTkbles;
 874             break;
 875           }
 876             if (jev >= fNpart || fNpart == -1) {
 877                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 878
 879                 fQ  += fPythia->GetVINT(51);
 880                 fX1 += fPythia->GetVINT(41);
 881                 fX2 += fPythia->GetVINT(42);
 882                 fTrialsRun += fTrials;
 883                 fNev++;
 884                 MakeHeader();
 885                 break;
 886             }
 887         }
 888     } // event loop
 889     SetHighWaterMark(nt);
 890 //  adjust weight due to kinematic selection
 891 //    AdjustWeights();
 892 //  get cross-section
 893     fXsection=fPythia->GetPARI(1);
 894 }
 895
 896 Int_t  AliGenPythia::GenerateMB()
 897 {
 898 //
 899 // Min Bias selection and other global selections
 900 //
 901     Int_t i, kf, nt, iparent;
 902     Int_t nc = 0;
 903     Double_t p[4];
 904     Double_t polar[3]   =   {0,0,0};
 905     Double_t origin[3]  =   {0,0,0};
 906 //  converts from mm/c to s
 907     const Float_t kconv=0.001/2.999792458e8;
 908
 909
 910
 911     Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
 912
 913
 914
 915     Int_t* pParent = new Int_t[np];
 916     for (i=0; i< np; i++) pParent[i] = -1;
 917      if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) {
 918         TParticle* jet1 = (TParticle *) fParticles.At(6);
 919         TParticle* jet2 = (TParticle *) fParticles.At(7);
 920         if (!CheckTrigger(jet1, jet2)) {
 921           delete [] pParent;
 922           return 0;
 923         }
 924     }
 925
 926     // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
 927     if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
 928
 929       Bool_t ok = kFALSE;
 930
 931       Int_t pdg  = 0;
 932       if (fFragPhotonInCalo) pdg = 22   ; // Photon
 933       else if (fPi0InCalo)   pdg = 111 ;    // Pi0
 934
 935       for (i=0; i< np; i++) {
 936         TParticle* iparticle = (TParticle *) fParticles.At(i);
 937         if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg &&
 938            iparticle->Pt() > fFragPhotonOrPi0MinPt){
 939           Int_t imother = iparticle->GetFirstMother() - 1;
 940           TParticle* pmother = (TParticle *) fParticles.At(imother);
 941           if(pdg == 111 ||
 942              (pdg == 22 && pmother->GetStatusCode() != 11)) //No photon from hadron decay
 943             {
 944               Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
 945               Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
 946               if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
 947                  (fCheckPHOS    && IsInPHOS(phi,eta)) )
 948                 ok =kTRUE;
 949             }
 950         }
 951       }
 952       if(!ok) {
 953           delete[] pParent;
 954           return 0;
 955       }
 956     }
 957
 958     // Select beauty jets with electron in EMCAL
 959     if (fProcess == kPyBeautyJets && fEleInEMCAL) {
 960
 961       Bool_t ok = kFALSE;
 962
 963       Int_t pdg  = 11; //electron
 964
 965       Float_t pt  = 0.;
 966       Float_t eta = 0.;
 967       Float_t phi = 0.;
 968       for (i=0; i< np; i++) {
 969         TParticle* iparticle = (TParticle *) fParticles.At(i);
 970         if(iparticle->GetStatusCode()==1 && TMath::Abs(iparticle->GetPdgCode())==pdg &&
 971            iparticle->Pt() > fElectronMinPt){
 972           pt = iparticle->Pt();
 973           phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
 974           eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
 975           if(IsInEMCAL(phi,eta))
 976             ok =kTRUE;
 977         }
 978       }
 979       if(!ok) {
 980         delete[] pParent;
 981         return 0;
 982       }
 983
 984       AliDebug(5,Form("Found an electron jet (pt,eta,phi) = (%f,%f,%f)",pt,eta,phi));
 985     }
 986     // Check for diffraction
 987     if(fkTuneForDiff) {
 988       if(fItune==320) {
 989         if(!CheckDiffraction()) {
 990           delete [] pParent;
 991           return 0;
 992         }
 993       }
 994     }
 995
 996     // Check for minimum multiplicity
 997     if (fTriggerMultiplicity > 0) {
 998       Int_t multiplicity = 0;
 999       for (i = 0; i < np; i++) {
1000         TParticle *  iparticle = (TParticle *) fParticles.At(i);
1001
1002         Int_t statusCode = iparticle->GetStatusCode();
1003
1004         // Initial state particle
1005         if (statusCode != 1)
1006           continue;
1007         // eta cut
1008         if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
1009           continue;
1010         // pt cut
1011         if (iparticle->Pt() < fTriggerMultiplicityPtMin)
1012             continue;
1013
1014         TParticlePDG* pdgPart = iparticle->GetPDG();
1015         if (pdgPart && pdgPart->Charge() == 0)
1016           continue;
1017
1018         ++multiplicity;
1019       }
1020
1021       if (multiplicity < fTriggerMultiplicity) {
1022         delete [] pParent;
1023         return 0;
1024       }
1025       Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
1026     }
1027
1028      // Select events with a photon  pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
1029     if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
1030
1031       Bool_t okd = kFALSE;
1032
1033       Int_t pdg  = 22;
1034       Int_t iphcand = -1;
1035       for (i=0; i< np; i++) {
1036          TParticle* iparticle = (TParticle *) fParticles.At(i);
1037          Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
1038          Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
1039
1040          if(iparticle->GetStatusCode() == 1
1041             && iparticle->GetPdgCode() == pdg
1042             && iparticle->Pt() > fPhotonMinPt
1043             && eta < fPHOSEta){
1044
1045             // first check if the photon is in PHOS phi
1046             if(IsInPHOS(phi,eta)){
1047                 okd = kTRUE;
1048                 break;
1049             }
1050             if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi
1051
1052          }
1053       }
1054
1055       if(!okd && iphcand != -1) // execute rotation in phi
1056           RotatePhi(iphcand,okd);
1057
1058       if(!okd) {
1059           delete [] pParent;
1060           return 0;
1061       }
1062     }
1063
1064     if (fTriggerParticle) {
1065         Bool_t triggered = kFALSE;
1066         for (i = 0; i < np; i++) {
1067             TParticle *  iparticle = (TParticle *) fParticles.At(i);
1068             kf = CheckPDGCode(iparticle->GetPdgCode());
1069             if (kf != fTriggerParticle) continue;
1070             if (iparticle->Pt() == 0.) continue;
1071             if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
1072             triggered = kTRUE;
1073             break;
1074         }
1075         if (!triggered) {
1076           delete [] pParent;
1077           return 0;
1078         }
1079     }
1080
1081
1082     // Check if there is a ccbar or bbbar pair with at least one of the two
1083     // in fYMin < y < fYMax
1084
1085     if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
1086       TParticle *partCheck;
1087       TParticle *mother;
1088       Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
1089       Bool_t  theChild=kFALSE;
1090       Bool_t  triggered=kFALSE;
1091       Float_t y;
1092       Int_t   pdg,mpdg,mpdgUpperFamily;
1093       for(i=0; i<np; i++) {
1094         partCheck = (TParticle*)fParticles.At(i);
1095         pdg = partCheck->GetPdgCode();
1096         if(TMath::Abs(pdg) == fFlavorSelect) { // quark
1097           if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1098           y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1099                              (partCheck->Energy()-partCheck->Pz()+1.e-13));
1100           if(y>fYMin && y<fYMax) inYcut=kTRUE;
1101         }
1102         if(fTriggerParticle) {
1103           if(TMath::Abs(pdg)==fTriggerParticle) triggered=kTRUE;
1104         }
1105         if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1106           Int_t mi = partCheck->GetFirstMother() - 1;
1107           if(mi<0) continue;
1108           mother = (TParticle*)fParticles.At(mi);
1109           mpdg=TMath::Abs(mother->GetPdgCode());
1110           mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
1111           if ( ParentSelected(mpdg) ||
1112               (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1113             if (KinematicSelection(partCheck,1)) {
1114               theChild=kTRUE;
1115             }
1116           }
1117         }
1118       }
1119       if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
1120         delete[] pParent;
1121         return 0;
1122       }
1123       if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1124         delete[] pParent;
1125         return 0;
1126       }
1127       if(fTriggerParticle && !triggered) { // particle requested is not produced
1128         delete[] pParent;
1129         return 0;
1130       }
1131
1132     }
1133
1134     //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
1135     if ( (fProcess == kPyW ||
1136           fProcess == kPyZ ||
1137           fProcess == kPyMbDefault ||
1138           fProcess == kPyMb ||
1139           fProcess == kPyMbAtlasTuneMC09 ||
1140           fProcess == kPyMbWithDirectPhoton ||
1141           fProcess == kPyMbNonDiffr)
1142          && (fCutOnChild == 1) ) {
1143       if ( !CheckKinematicsOnChild() ) {
1144         delete[] pParent;
1145         return 0;
1146       }
1147     }
1148
1149
1150     for (i = 0; i < np; i++) {
1151         Int_t trackIt = 0;
1152         TParticle *  iparticle = (TParticle *) fParticles.At(i);
1153         kf = CheckPDGCode(iparticle->GetPdgCode());
1154         Int_t ks = iparticle->GetStatusCode();
1155         Int_t km = iparticle->GetFirstMother();
1156         if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
1157             (ks != 1) ||
1158             ((fProcess == kPyJets || fProcess == kPyBeautyJets) && ks == 21 && km == 0 && i>1)) {
1159             nc++;
1160             if (ks == 1) trackIt = 1;
1161             Int_t ipa = iparticle->GetFirstMother()-1;
1162
1163             iparent = (ipa > -1) ? pParent[ipa] : -1;
1164
1165 //
1166 // store track information
1167             p[0] = iparticle->Px();
1168             p[1] = iparticle->Py();
1169             p[2] = iparticle->Pz();
1170             p[3] = iparticle->Energy();
1171
1172
1173             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1174             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1175             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1176
1177             Float_t tof = fEventTime + kconv * iparticle->T();
1178
1179             PushTrack(fTrackIt*trackIt, iparent, kf,
1180                       p[0], p[1], p[2], p[3],
1181                       origin[0], origin[1], origin[2], tof,
1182                       polar[0], polar[1], polar[2],
1183                       kPPrimary, nt, 1., ks);
1184             fNprimaries++;
1185             KeepTrack(nt);
1186             pParent[i] = nt;
1187             SetHighWaterMark(nt);
1188
1189         } // select particle
1190     } // particle loop
1191
1192     delete[] pParent;
1193
1194     return 1;
1195 }
1196
1197
1198 void AliGenPythia::FinishRun()
1199 {
1200 // Print x-section summary
1201     fPythia->Pystat(1);
1202
1203     if (fNev > 0.) {
1204         fQ  /= fNev;
1205         fX1 /= fNev;
1206         fX2 /= fNev;
1207     }
1208
1209     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1210     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1211 }
1212
1213 void AliGenPythia::AdjustWeights() const
1214 {
1215 // Adjust the weights after generation of all events
1216 //
1217     if (gAlice) {
1218         TParticle *part;
1219         Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1220         for (Int_t i=0; i<ntrack; i++) {
1221             part= gAlice->GetMCApp()->Particle(i);
1222             part->SetWeight(part->GetWeight()*fKineBias);
1223         }
1224     }
1225 }
1226
1227 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
1228 {
1229 // Treat protons as inside nuclei with mass numbers a1 and a2
1230
1231     fAProjectile = a1;
1232     fATarget     = a2;
1233     fNucPdf      = pdfset;  // 0 EKS98 1 EPS08
1234     fSetNuclei   = kTRUE;
1235 }
1236
1237
1238 void AliGenPythia::MakeHeader()
1239 {
1240 //
1241 // Make header for the simulated event
1242 //
1243   if (gAlice) {
1244     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1245         gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1246   }
1247
1248 // Builds the event header, to be called after each event
1249     if (fHeader) delete fHeader;
1250     fHeader = new AliGenPythiaEventHeader("Pythia");
1251 //
1252 // Event type
1253     if(fProcDiff>0){
1254       //      if(fProcDiff == 92 || fProcDiff == 93) printf("\n\n\n\n\n");
1255       //      printf("fPythia->GetMSTI(1) = %d   fProcDiff = %d\n",fPythia->GetMSTI(1), fProcDiff);
1256     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fProcDiff);
1257     }
1258     else
1259     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1260 //
1261 // Number of trials
1262     ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1263 //
1264 // Event Vertex
1265     fHeader->SetPrimaryVertex(fVertex);
1266     fHeader->SetInteractionTime(fEventTime);
1267 //
1268 // Number of primaries
1269     fHeader->SetNProduced(fNprimaries);
1270 //
1271 // Jets that have triggered
1272
1273     //Need to store jets for b-jet studies too!
1274     if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi)
1275     {
1276         Int_t ntrig, njet;
1277         Float_t jets[4][10];
1278         GetJets(njet, ntrig, jets);
1279
1280
1281         for (Int_t i = 0; i < ntrig; i++) {
1282             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1283                                                         jets[3][i]);
1284         }
1285     }
1286 //
1287 // Copy relevant information from external header, if present.
1288 //
1289     Float_t uqJet[4];
1290
1291     if (fRL) {
1292         AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1293         for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1294         {
1295             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
1296
1297
1298             exHeader->TriggerJet(i, uqJet);
1299             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1300         }
1301     }
1302 //
1303 // Store quenching parameters
1304 //
1305     if (fQuench){
1306         Double_t z[4] = {0.};
1307         Double_t xp = 0.;
1308         Double_t yp = 0.;
1309
1310         if (fQuench == 1) {
1311             // Pythia::Quench()
1312             fPythia->GetQuenchingParameters(xp, yp, z);
1313         } else if (fQuench == 2){
1314             // Pyquen
1315             Double_t r1 = PARIMP.rb1;
1316             Double_t r2 = PARIMP.rb2;
1317             Double_t b  = PARIMP.b1;
1318             Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1319             Double_t phi = PARIMP.psib1;
1320             xp = r * TMath::Cos(phi);
1321             yp = r * TMath::Sin(phi);
1322
1323         } else if (fQuench == 4) {
1324             // QPythia
1325             Double_t xy[2];
1326             Double_t i0i1[2];
1327             AliFastGlauber::Instance()->GetSavedXY(xy);
1328             AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
1329             xp = xy[0];
1330             yp = xy[1];
1331             ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
1332         }
1333
1334             ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1335             ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1336     }
1337 //
1338 // Store pt^hard
1339     ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1340 //
1341 //  Pass header
1342 //
1343     AddHeader(fHeader);
1344     fHeader = 0x0;
1345 }
1346
1347 Bool_t AliGenPythia::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
1348 {
1349 // Check the kinematic trigger condition
1350 //
1351     Double_t eta[2];
1352     eta[0] = jet1->Eta();
1353     eta[1] = jet2->Eta();
1354     Double_t phi[2];
1355     phi[0] = jet1->Phi();
1356     phi[1] = jet2->Phi();
1357     Int_t    pdg[2];
1358     pdg[0] = jet1->GetPdgCode();
1359     pdg[1] = jet2->GetPdgCode();
1360     Bool_t   triggered = kFALSE;
1361
1362     if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess ==  kPyBeautyppMNRwmi) {
1363         Int_t njets = 0;
1364         Int_t ntrig = 0;
1365         Float_t jets[4][10];
1366 //
1367 // Use Pythia clustering on parton level to determine jet axis
1368 //
1369         GetJets(njets, ntrig, jets);
1370
1371         if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1372 //
1373     } else {
1374         Int_t ij = 0;
1375         Int_t ig = 1;
1376         if (pdg[0] == kGamma) {
1377             ij = 1;
1378             ig = 0;
1379         }
1380         //Check eta range first...
1381         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
1382             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1383         {
1384             //Eta is okay, now check phi range
1385             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
1386                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1387             {
1388                 triggered = kTRUE;
1389             }
1390         }
1391     }
1392     return triggered;
1393 }
1394
1395
1396
1397 Bool_t AliGenPythia::CheckKinematicsOnChild(){
1398 //
1399 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1400 //
1401     Bool_t checking = kFALSE;
1402     Int_t j, kcode, ks, km;
1403     Int_t nPartAcc = 0; //number of particles in the acceptance range
1404     Int_t numberOfAcceptedParticles = 1;
1405     if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1406     Int_t npart = fParticles.GetEntriesFast();
1407
1408     for (j = 0; j<npart; j++) {
1409         TParticle *  jparticle = (TParticle *) fParticles.At(j);
1410         kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1411         ks = jparticle->GetStatusCode();
1412         km = jparticle->GetFirstMother();
1413
1414         if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
1415             nPartAcc++;
1416         }
1417         if( numberOfAcceptedParticles <= nPartAcc){
1418           checking = kTRUE;
1419           break;
1420         }
1421     }
1422
1423     return checking;
1424 }
1425
1426 void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1427 {
1428   //
1429   // Load event into Pythia Common Block
1430   //
1431
1432   Int_t npart = stack -> GetNprimary();
1433   Int_t n0 = 0;
1434
1435   if (!flag) {
1436     (fPythia->GetPyjets())->N = npart;
1437   } else {
1438     n0 = (fPythia->GetPyjets())->N;
1439     (fPythia->GetPyjets())->N = n0 + npart;
1440   }
1441
1442
1443   for (Int_t part = 0; part < npart; part++) {
1444     TParticle *mPart = stack->Particle(part);
1445
1446     Int_t kf     =  mPart->GetPdgCode();
1447     Int_t ks     =  mPart->GetStatusCode();
1448     Int_t idf    =  mPart->GetFirstDaughter();
1449     Int_t idl    =  mPart->GetLastDaughter();
1450
1451     if (reHadr) {
1452             if (ks == 11 || ks == 12) {
1453         ks  -= 10;
1454         idf  = -1;
1455         idl  = -1;
1456             }
1457     }
1458
1459     Float_t px = mPart->Px();
1460     Float_t py = mPart->Py();
1461     Float_t pz = mPart->Pz();
1462     Float_t e  = mPart->Energy();
1463     Float_t m  = mPart->GetCalcMass();
1464
1465
1466     (fPythia->GetPyjets())->P[0][part+n0] = px;
1467     (fPythia->GetPyjets())->P[1][part+n0] = py;
1468     (fPythia->GetPyjets())->P[2][part+n0] = pz;
1469     (fPythia->GetPyjets())->P[3][part+n0] = e;
1470     (fPythia->GetPyjets())->P[4][part+n0] = m;
1471
1472     (fPythia->GetPyjets())->K[1][part+n0] = kf;
1473     (fPythia->GetPyjets())->K[0][part+n0] = ks;
1474     (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1475     (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1476     (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1477   }
1478 }
1479
1480 void  AliGenPythia::LoadEvent(const TObjArray* stack, Int_t flag, Int_t reHadr)
1481 {
1482   //
1483   // Load event into Pythia Common Block
1484   //
1485
1486   Int_t npart = stack -> GetEntries();
1487   Int_t n0 = 0;
1488
1489   if (!flag) {
1490     (fPythia->GetPyjets())->N = npart;
1491   } else {
1492     n0 = (fPythia->GetPyjets())->N;
1493     (fPythia->GetPyjets())->N = n0 + npart;
1494   }
1495
1496
1497   for (Int_t part = 0; part < npart; part++) {
1498     TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
1499     if (!mPart) continue;
1500
1501     Int_t kf     =  mPart->GetPdgCode();
1502     Int_t ks     =  mPart->GetStatusCode();
1503     Int_t idf    =  mPart->GetFirstDaughter();
1504     Int_t idl    =  mPart->GetLastDaughter();
1505
1506     if (reHadr) {
1507         if (ks == 11 || ks == 12) {
1508             ks  -= 10;
1509             idf  = -1;
1510             idl  = -1;
1511         }
1512     }
1513
1514     Float_t px = mPart->Px();
1515     Float_t py = mPart->Py();
1516     Float_t pz = mPart->Pz();
1517     Float_t e  = mPart->Energy();
1518     Float_t m  = mPart->GetCalcMass();
1519
1520
1521     (fPythia->GetPyjets())->P[0][part+n0] = px;
1522     (fPythia->GetPyjets())->P[1][part+n0] = py;
1523     (fPythia->GetPyjets())->P[2][part+n0] = pz;
1524     (fPythia->GetPyjets())->P[3][part+n0] = e;
1525     (fPythia->GetPyjets())->P[4][part+n0] = m;
1526
1527     (fPythia->GetPyjets())->K[1][part+n0] = kf;
1528     (fPythia->GetPyjets())->K[0][part+n0] = ks;
1529     (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1530     (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1531     (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1532   }
1533 }
1534
1535
1536 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1537 {
1538 //
1539 //  Calls the Pythia jet finding algorithm to find jets in the current event
1540 //
1541 //
1542 //
1543 //  Save jets
1544     Int_t n     = fPythia->GetN();
1545
1546 //
1547 //  Run Jet Finder
1548     fPythia->Pycell(njets);
1549     Int_t i;
1550     for (i = 0; i < njets; i++) {
1551         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1552         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1553         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1554         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1555
1556         jets[0][i] = px;
1557         jets[1][i] = py;
1558         jets[2][i] = pz;
1559         jets[3][i] = e;
1560     }
1561 }
1562
1563
1564
1565 void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1566 {
1567 //
1568 //  Calls the Pythia clustering algorithm to find jets in the current event
1569 //
1570     Int_t n     = fPythia->GetN();
1571     nJets       = 0;
1572     nJetsTrig   = 0;
1573     if (fJetReconstruction == kCluster) {
1574 //
1575 //  Configure cluster algorithm
1576 //
1577         fPythia->SetPARU(43, 2.);
1578         fPythia->SetMSTU(41, 1);
1579 //
1580 //  Call cluster algorithm
1581 //
1582         fPythia->Pyclus(nJets);
1583 //
1584 //  Loading jets from common block
1585 //
1586     } else {
1587
1588 //
1589 //  Run Jet Finder
1590         fPythia->Pycell(nJets);
1591     }
1592
1593     Int_t i;
1594     for (i = 0; i < nJets; i++) {
1595         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1596         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1597         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1598         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1599         Float_t pt    = TMath::Sqrt(px * px + py * py);
1600         Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);
1601         Float_t theta = TMath::ATan2(pt,pz);
1602         Float_t et    = e * TMath::Sin(theta);
1603         Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
1604         if (
1605             eta > fEtaMinJet && eta < fEtaMaxJet &&
1606             phi > fPhiMinJet && phi < fPhiMaxJet &&
1607             et  > fEtMinJet  && et  < fEtMaxJet
1608             )
1609         {
1610             jets[0][nJetsTrig] = px;
1611             jets[1][nJetsTrig] = py;
1612             jets[2][nJetsTrig] = pz;
1613             jets[3][nJetsTrig] = e;
1614             nJetsTrig++;
1615 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1616         } else {
1617 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1618         }
1619     }
1620 }
1621
1622 void AliGenPythia::GetSubEventTime()
1623 {
1624   // Calculates time of the next subevent
1625   fEventTime = 0.;
1626   if (fEventsTime) {
1627     TArrayF &array = *fEventsTime;
1628     fEventTime = array[fCurSubEvent++];
1629   }
1630   //  printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1631   return;
1632 }
1633
1634 Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta) const
1635 {
1636   // Is particle in EMCAL acceptance?
1637   // phi in degrees, etamin=-etamax
1638   if(phi > fEMCALMinPhi  && phi < fEMCALMaxPhi &&
1639      eta < fEMCALEta  )
1640     return kTRUE;
1641   else
1642     return kFALSE;
1643 }
1644
1645 Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta) const
1646 {
1647   // Is particle in PHOS acceptance?
1648   // Acceptance slightly larger considered.
1649   // phi in degrees, etamin=-etamax
1650   if(phi > fPHOSMinPhi  && phi < fPHOSMaxPhi &&
1651      eta < fPHOSEta  )
1652     return kTRUE;
1653   else
1654     return kFALSE;
1655 }
1656
1657 void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd)
1658 {
1659   //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1660   Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1661   Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1662   Double_t phiPHOS = (AliPythiaRndm::GetPythiaRandom())->Uniform(phiPHOSmin,phiPHOSmax);
1663
1664   //calculate deltaphi
1665   TParticle* ph = (TParticle *) fParticles.At(iphcand);
1666   Double_t phphi = ph->Phi();
1667   Double_t deltaphi = phiPHOS - phphi;
1668
1669
1670
1671   //loop for all particles and produce the phi rotation
1672   Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1673   Double_t oldphi, newphi;
1674   Double_t newVx, newVy, r, vZ, time;
1675   Double_t newPx, newPy, pt, pz, e;
1676   for(Int_t i=0; i< np; i++) {
1677       TParticle* iparticle = (TParticle *) fParticles.At(i);
1678       oldphi = iparticle->Phi();
1679       newphi = oldphi + deltaphi;
1680       if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1681       if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1682
1683       r = iparticle->R();
1684       newVx = r * TMath::Cos(newphi);
1685       newVy = r * TMath::Sin(newphi);
1686       vZ   = iparticle->Vz(); // don't transform
1687       time = iparticle->T(); // don't transform
1688
1689       pt = iparticle->Pt();
1690       newPx = pt * TMath::Cos(newphi);
1691       newPy = pt * TMath::Sin(newphi);
1692       pz = iparticle->Pz(); // don't transform
1693       e  = iparticle->Energy(); // don't transform
1694
1695       // apply rotation
1696       iparticle->SetProductionVertex(newVx, newVy, vZ, time);
1697       iparticle->SetMomentum(newPx, newPy, pz, e);
1698
1699   } //end particle loop
1700
1701    // now let's check that we put correctly the candidate photon in PHOS
1702    Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1703    Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1704    if(IsInPHOS(phi,eta))
1705       okdd = kTRUE;
1706 }
1707
1708
1709
1710 Bool_t AliGenPythia::CheckDiffraction()
1711 {
1712   // use this method only with Perugia-0 tune!
1713
1714   //  printf("AAA\n");
1715
1716    Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1717
1718    Int_t iPart1=-1;
1719    Int_t iPart2=-1;
1720
1721    Double_t y1 = 1e10;
1722    Double_t y2 = -1e10;
1723
1724   const Int_t kNstable=20;
1725   const Int_t pdgStable[20] = {
1726     22,             // Photon
1727     11,             // Electron
1728     12,             // Electron Neutrino
1729     13,             // Muon
1730     14,             // Muon Neutrino
1731     15,             // Tau
1732     16,             // Tau Neutrino
1733     211,            // Pion
1734     321,            // Kaon
1735     311,            // K0
1736     130,            // K0s
1737     310,            // K0l
1738     2212,           // Proton
1739     2112,           // Neutron
1740     3122,           // Lambda_0
1741     3112,           // Sigma Minus
1742     3222,           // Sigma Plus
1743     3312,           // Xsi Minus
1744     3322,           // Xsi0
1745     3334            // Omega
1746   };
1747
1748      for (Int_t i = 0; i < np; i++) {
1749         TParticle *  part = (TParticle *) fParticles.At(i);
1750
1751         Int_t statusCode = part->GetStatusCode();
1752
1753         // Initial state particle
1754         if (statusCode != 1)
1755           continue;
1756
1757         Int_t pdg = TMath::Abs(part->GetPdgCode());
1758         Bool_t isStable = kFALSE;
1759         for (Int_t i1 = 0; i1 < kNstable; i1++) {
1760           if (pdg == pdgStable[i1]) {
1761             isStable = kTRUE;
1762             break;
1763           }
1764         }
1765         if(!isStable)
1766           continue;
1767
1768         Double_t y = part->Y();
1769
1770         if (y < y1)
1771           {
1772             y1 = y;
1773             iPart1 = i;
1774           }
1775         if (y > y2)
1776         {
1777           y2 = y;
1778           iPart2 = i;
1779         }
1780      }
1781
1782      if(iPart1<0 || iPart2<0) return kFALSE;
1783
1784      y1=TMath::Abs(y1);
1785      y2=TMath::Abs(y2);
1786
1787      TParticle *  part1 = (TParticle *) fParticles.At(iPart1);
1788      TParticle *  part2 = (TParticle *) fParticles.At(iPart2);
1789
1790      Int_t pdg1 = part1->GetPdgCode();
1791      Int_t pdg2 = part2->GetPdgCode();
1792
1793
1794      Int_t iPart = -1;
1795      if (pdg1 == 2212 && pdg2 == 2212)
1796        {
1797          if(y1 > y2)
1798            iPart = iPart1;
1799          else if(y1 < y2)
1800            iPart = iPart2;
1801          else {
1802            iPart = iPart1;
1803            if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)>0.5) iPart = iPart2;
1804          }
1805        }
1806      else if (pdg1 == 2212)
1807        iPart = iPart1;
1808      else if (pdg2 == 2212)
1809        iPart = iPart2;
1810
1811
1812
1813
1814
1815      Double_t M=-1.;
1816      if(iPart>0) {
1817        TParticle *  part = (TParticle *) fParticles.At(iPart);
1818        Double_t E= part->Energy();
1819        Double_t P= part->P();
1820        M= TMath::Sqrt((fEnergyCMS-E-P)*(fEnergyCMS-E+P));
1821      }
1822
1823 const Int_t nbin=120;
1824 Double_t bin[]={
1825 1.080000, 1.274258, 1.468516, 1.662773, 1.857031, 2.051289,
1826 2.245547, 2.439805, 2.634062, 2.828320, 3.022578, 3.216836,
1827 3.411094, 3.605352, 3.799609, 3.993867, 4.188125, 4.382383,
1828 4.576641, 4.770898, 4.965156, 5.547930, 6.130703, 6.713477,
1829 7.296250, 7.879023, 8.461797, 9.044570, 9.627344, 10.210117,
1830 10.792891, 11.375664, 11.958437, 12.541211, 13.123984, 13.706758,
1831 14.289531, 14.872305, 15.455078, 16.037852, 16.620625, 17.203398,
1832 17.786172, 18.368945, 18.951719, 19.534492, 20.117266, 20.700039,
1833 21.282812, 21.865586, 22.448359, 23.031133, 23.613906, 24.196680,
1834 24.779453, 25.362227, 25.945000, 26.527773, 27.110547, 27.693320,
1835 28.276094, 28.858867, 29.441641, 30.024414, 30.607187, 31.189961,
1836 31.772734, 32.355508, 32.938281, 33.521055, 34.103828, 34.686602,
1837 35.269375, 35.852148, 36.434922, 37.017695, 37.600469, 38.183242,
1838 38.766016, 39.348789, 39.931562, 40.514336, 41.097109, 41.679883,
1839 42.262656, 42.845430, 43.428203, 44.010977, 44.593750, 45.176523,
1840 45.759297, 46.342070, 46.924844, 47.507617, 48.090391, 48.673164,
1841 49.255937, 49.838711, 50.421484, 57.220508, 64.019531, 70.818555,
1842 77.617578, 84.416602, 91.215625, 98.014648, 104.813672, 111.612695,
1843 118.411719, 125.210742, 132.009766, 138.808789, 145.607812, 152.406836,
1844 159.205859, 166.004883, 172.803906, 179.602930, 186.401953, 193.200977,
1845 200.000000};
1846 Double_t w[]={
1847 1.000000, 0.275812, 0.267308, 0.263268, 0.263828, 0.263039,
1848 0.260546, 0.259344, 0.255477, 0.253782, 0.253562, 0.252492,
1849 0.251076, 0.247862, 0.248933, 0.243599, 0.244255, 0.238506,
1850 0.239546, 0.237845, 0.237977, 0.229379, 0.224771, 0.217581,
1851 0.208860, 0.207241, 0.198955, 0.196449, 0.193827, 0.190220,
1852 0.184715, 0.183067, 0.178325, 0.175887, 0.171299, 0.168718,
1853 0.167453, 0.165153, 0.163457, 0.159637, 0.156855, 0.155488,
1854 0.154545, 0.155968, 0.150488, 0.148797, 0.145358, 0.146196,
1855 0.145891, 0.142752, 0.145072, 0.141265, 0.141857, 0.138462,
1856 0.142992, 0.141357, 0.139391, 0.139629, 0.135197, 0.135731,
1857 0.133194, 0.137190, 0.135745, 0.134522, 0.136094, 0.130405,
1858 0.128371, 0.131680, 0.130591, 0.133539, 0.129370, 0.128254,
1859 0.128262, 0.131088, 0.128294, 0.130070, 0.124553, 0.131489,
1860 0.128038, 0.122997, 0.130699, 0.125630, 0.124746, 0.123679,
1861 0.127864, 0.125776, 0.126272, 0.121492, 0.124929, 0.122221,
1862 0.127017, 0.123706, 0.122701, 0.123731, 0.122219, 0.120783,
1863 0.120324, 0.120228, 0.123615, 0.120589, 0.119549, 0.118836,
1864 0.118146, 0.120175, 0.122031, 0.122076, 0.122849, 0.122627,
1865 0.126281, 0.127696, 0.129084, 0.130001, 0.134062, 0.134786,
1866 0.137286, 0.136625, 0.139091, 0.143692, 0.144781, 0.146768};
1867
1868  Double_t wSD=1.;
1869  Double_t wDD=0.178783;
1870  //Double_t wND=0.220200;
1871  Double_t wND=0.220200+0.001;
1872
1873  if(M>-1 && M<bin[0]) return kFALSE;
1874  if(M>bin[nbin]) M=-1;
1875
1876  Int_t procType=fPythia->GetMSTI(1);
1877  Int_t proc0=2;
1878  if(procType== 94) proc0=1;
1879  if(procType== 92 || procType== 93) proc0=0;
1880
1881
1882  // printf("M = %f   bin[nbin] = %f\n",M, bin[nbin]);
1883
1884  Int_t proc=2;
1885  if(M>0) proc=0;
1886  else if(proc0==1) proc=1;
1887
1888  if(proc==0 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wSD) return kFALSE;
1889  if(proc==1 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wDD) return kFALSE;
1890  if(proc==2 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wND) return kFALSE;
1891
1892
1893     //     if(proc==1 || proc==2) return kFALSE;
1894
1895     if(proc!=0) {
1896       if(proc0!=0) fProcDiff = procType;
1897       else       fProcDiff = 95;
1898       return kTRUE;
1899     }
1900
1901     Int_t ibin=nbin-1;
1902     for(Int_t i=1; i<=nbin; i++)
1903       if(M<=bin[i]) {
1904         ibin=i-1;
1905         //      printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
1906         break;
1907       }
1908
1909     //    printf("w[ibin] = %f\n", w[ibin]);
1910
1911     if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)> w[ibin]) return kFALSE;
1912
1913     //    printf("iPart = %d\n", iPart);
1914
1915     if(iPart==iPart1) fProcDiff=93;
1916     else if(iPart==iPart2) fProcDiff=92;
1917     else {
1918       printf("EROOR:  iPart!=iPart1 && iPart!=iPart2\n");
1919
1920     }
1921
1922     return kTRUE;
1923 }