PYTHIA6/AliPythia6/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
  50 ClassImp(AliGenPythia)
  51
  52
  53 AliGenPythia::AliGenPythia():
  54     AliGenMC(),
  55     fProcess(kPyCharm),
  56     fItune(-1),
  57     fStrucFunc(kCTEQ5L),
  58     fKineBias(0.),
  59     fTrials(0),
  60     fTrialsRun(0),
  61     fQ(0.),
  62     fX1(0.),
  63     fX2(0.),
  64     fEventTime(0.),
  65     fInteractionRate(0.),
  66     fTimeWindow(0.),
  67     fCurSubEvent(0),
  68     fEventsTime(0),
  69     fNev(0),
  70     fFlavorSelect(0),
  71     fXsection(0.),
  72     fPythia(0),
  73     fWeightPower(0.),
  74     fPtHardMin(0.),
  75     fPtHardMax(1.e4),
  76     fYHardMin(-1.e10),
  77     fYHardMax(1.e10),
  78     fGinit(1),
  79     fGfinal(1),
  80     fCRoff(0),
  81     fHadronisation(1),
  82     fPatchOmegaDalitz(0),
  83     fDecayerExodus(0),
  84     fNpartons(0),
  85     fReadFromFile(0),
  86     fReadLHEF(0),
  87     fQuench(0),
  88     fQhat(0.),
  89     fLength(0.),
  90     fpyquenT(1.),
  91     fpyquenTau(0.1),
  92     fpyquenNf(0),
  93     fpyquenEloss(0),
  94     fpyquenAngle(0),
  95     fImpact(0.),
  96     fPtKick(1.),
  97     fFullEvent(kTRUE),
  98     fDecayer(new AliDecayerPythia()),
  99     fDebugEventFirst(-1),
 100     fDebugEventLast(-1),
 101     fEtMinJet(0.),
 102     fEtMaxJet(1.e4),
 103     fEtaMinJet(-20.),
 104     fEtaMaxJet(20.),
 105     fPhiMinJet(0.),
 106     fPhiMaxJet(2.* TMath::Pi()),
 107     fJetReconstruction(kCell),
 108     fEtaMinGamma(-20.),
 109     fEtaMaxGamma(20.),
 110     fPhiMinGamma(0.),
 111     fPhiMaxGamma(2. * TMath::Pi()),
 112     fUseYCutHQ(kFALSE),
 113     fYMinHQ(-20.),
 114     fYMaxHQ(20.),
 115     fPycellEtaMax(2.),
 116     fPycellNEta(274),
 117     fPycellNPhi(432),
 118     fPycellThreshold(0.),
 119     fPycellEtSeed(4.),
 120     fPycellMinEtJet(10.),
 121     fPycellMaxRadius(1.),
 122     fStackFillOpt(kFlavorSelection),
 123     fFeedDownOpt(kTRUE),
 124     fFragmentation(kTRUE),
 125     fSetNuclei(kFALSE),
 126     fUseNuclearPDF(kFALSE),
 127     fUseLorentzBoost(kTRUE),
 128     fNewMIS(kFALSE),
 129     fHFoff(kFALSE),
 130     fNucPdf(0),
 131     fTriggerParticle(0),
 132     fTriggerEta(0.9),
 133     fTriggerY(999.),
 134     fTriggerEtaMin(0.9),
 135     fTriggerMinPt(-1),
 136     fTriggerMaxPt(1000),
 137     fTriggerMultiplicity(0),
 138     fTriggerMultiplicityEta(0),
 139     fTriggerMultiplicityEtaMin(0),
 140     fTriggerMultiplicityEtaMax(0),
 141     fTriggerMultiplicityPtMin(0),
 142     fCountMode(kCountAll),
 143     fHeader(0),
 144     fRL(0),
 145     fkFileName(0),
 146     fkNameLHEF(0),
 147     fFragPhotonInCalo(kFALSE),
 148     fHadronInCalo(kFALSE) ,
 149     fPi0InCalo(kFALSE) ,
 150     fEtaInCalo(kFALSE) ,
 151     fPhotonInCalo(kFALSE), // not in use
 152     fDecayPhotonInCalo(kFALSE),
 153     fForceNeutralMeson2PhotonDecay(kFALSE),
 154     fEleInCalo(kFALSE),
 155     fEleInEMCAL(kFALSE), // not in use
 156     fCheckBarrel(kFALSE),
 157     fCheckEMCAL(kFALSE),
 158     fCheckPHOS(kFALSE),
 159     fCheckPHOSeta(kFALSE),
 160     fPHOSRotateCandidate(-1),
 161     fTriggerParticleMinPt(0),
 162     fPhotonMinPt(0), // not in use
 163     fElectronMinPt(0), // not in use
 164     fPHOSMinPhi(219.),
 165     fPHOSMaxPhi(321.),
 166     fPHOSEta(0.13),
 167     fEMCALMinPhi(79.),
 168     fEMCALMaxPhi(191.),
 169     fEMCALEta(0.71),
 170     fkTuneForDiff(0),
 171     fProcDiff(0)
 172 {
 173 // Default Constructor
 174   fEnergyCMS = 5500.;
 175   if (!AliPythiaRndm::GetPythiaRandom())
 176       AliPythiaRndm::SetPythiaRandom(GetRandom());
 177 }
 178
 179 AliGenPythia::AliGenPythia(Int_t npart)
 180     :AliGenMC(npart),
 181      fProcess(kPyCharm),
 182      fItune(-1),
 183      fStrucFunc(kCTEQ5L),
 184      fKineBias(0.),
 185      fTrials(0),
 186      fTrialsRun(0),
 187      fQ(0.),
 188      fX1(0.),
 189      fX2(0.),
 190      fEventTime(0.),
 191      fInteractionRate(0.),
 192      fTimeWindow(0.),
 193      fCurSubEvent(0),
 194      fEventsTime(0),
 195      fNev(0),
 196      fFlavorSelect(0),
 197      fXsection(0.),
 198      fPythia(0),
 199      fWeightPower(0.),
 200      fPtHardMin(0.),
 201      fPtHardMax(1.e4),
 202      fYHardMin(-1.e10),
 203      fYHardMax(1.e10),
 204      fGinit(kTRUE),
 205      fGfinal(kTRUE),
 206      fCRoff(kFALSE),
 207      fHadronisation(kTRUE),
 208      fPatchOmegaDalitz(0),
 209      fDecayerExodus(0),
 210      fNpartons(0),
 211      fReadFromFile(kFALSE),
 212      fReadLHEF(0),
 213      fQuench(kFALSE),
 214      fQhat(0.),
 215      fLength(0.),
 216      fpyquenT(1.),
 217      fpyquenTau(0.1),
 218      fpyquenNf(0),
 219      fpyquenEloss(0),
 220      fpyquenAngle(0),
 221      fImpact(0.),
 222      fPtKick(1.),
 223      fFullEvent(kTRUE),
 224      fDecayer(new AliDecayerPythia()),
 225      fDebugEventFirst(-1),
 226      fDebugEventLast(-1),
 227      fEtMinJet(0.),
 228      fEtMaxJet(1.e4),
 229      fEtaMinJet(-20.),
 230      fEtaMaxJet(20.),
 231      fPhiMinJet(0.),
 232      fPhiMaxJet(2.* TMath::Pi()),
 233      fJetReconstruction(kCell),
 234      fEtaMinGamma(-20.),
 235      fEtaMaxGamma(20.),
 236      fPhiMinGamma(0.),
 237      fPhiMaxGamma(2. * TMath::Pi()),
 238      fUseYCutHQ(kFALSE),
 239      fYMinHQ(-20.),
 240      fYMaxHQ(20.),
 241      fPycellEtaMax(2.),
 242      fPycellNEta(274),
 243      fPycellNPhi(432),
 244      fPycellThreshold(0.),
 245      fPycellEtSeed(4.),
 246      fPycellMinEtJet(10.),
 247      fPycellMaxRadius(1.),
 248      fStackFillOpt(kFlavorSelection),
 249      fFeedDownOpt(kTRUE),
 250      fFragmentation(kTRUE),
 251      fSetNuclei(kFALSE),
 252      fUseNuclearPDF(kFALSE),
 253      fUseLorentzBoost(kTRUE),
 254      fNewMIS(kFALSE),
 255      fHFoff(kFALSE),
 256      fNucPdf(0),
 257      fTriggerParticle(0),
 258      fTriggerEta(0.9),
 259      fTriggerY(999.),
 260      fTriggerEtaMin(0.9),
 261      fTriggerMinPt(-1),
 262      fTriggerMaxPt(1000),
 263      fTriggerMultiplicity(0),
 264      fTriggerMultiplicityEta(0),
 265      fTriggerMultiplicityEtaMin(0),
 266      fTriggerMultiplicityEtaMax(0),
 267      fTriggerMultiplicityPtMin(0),
 268      fCountMode(kCountAll),
 269      fHeader(0),
 270      fRL(0),
 271      fkFileName(0),
 272      fkNameLHEF(0),
 273      fFragPhotonInCalo(kFALSE),
 274      fHadronInCalo(kFALSE) ,
 275      fPi0InCalo(kFALSE) ,
 276      fEtaInCalo(kFALSE) ,
 277      fPhotonInCalo(kFALSE), // not in use
 278      fDecayPhotonInCalo(kFALSE),
 279      fForceNeutralMeson2PhotonDecay(kFALSE),
 280      fEleInCalo(kFALSE),
 281      fEleInEMCAL(kFALSE), // not in use
 282      fCheckBarrel(kFALSE),
 283      fCheckEMCAL(kFALSE),
 284      fCheckPHOS(kFALSE),
 285      fCheckPHOSeta(kFALSE),
 286      fPHOSRotateCandidate(-1),
 287      fTriggerParticleMinPt(0),
 288      fPhotonMinPt(0), // not in use
 289      fElectronMinPt(0), // not in use
 290      fPHOSMinPhi(219.),
 291      fPHOSMaxPhi(321.),
 292      fPHOSEta(0.13),
 293      fEMCALMinPhi(79.),
 294      fEMCALMaxPhi(191.),
 295      fEMCALEta(0.71),
 296      fkTuneForDiff(0),
 297      fProcDiff(0)
 298 {
 299 // default charm production at 5. 5 TeV
 300 // semimuonic decay
 301 // structure function GRVHO
 302 //
 303     fEnergyCMS = 5500.;
 304     fName = "Pythia";
 305     fTitle= "Particle Generator using PYTHIA";
 306     SetForceDecay();
 307     // Set random number generator
 308     if (!AliPythiaRndm::GetPythiaRandom())
 309       AliPythiaRndm::SetPythiaRandom(GetRandom());
 310  }
 311
 312 AliGenPythia::~AliGenPythia()
 313 {
 314 // Destructor
 315   if(fEventsTime) delete fEventsTime;
 316 }
 317
 318 void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
 319 {
 320 // Generate pileup using user specified rate
 321     fInteractionRate = rate;
 322     fTimeWindow = timewindow;
 323     GeneratePileup();
 324 }
 325
 326 void AliGenPythia::GeneratePileup()
 327 {
 328 // Generate sub events time for pileup
 329     fEventsTime = 0;
 330     if(fInteractionRate == 0.) {
 331       Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
 332       return;
 333     }
 334
 335     Int_t npart = NumberParticles();
 336     if(npart < 0) {
 337       Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
 338       return;
 339     }
 340
 341     if(fEventsTime) delete fEventsTime;
 342     fEventsTime = new TArrayF(npart);
 343     TArrayF &array = *fEventsTime;
 344     for(Int_t ipart = 0; ipart < npart; ipart++)
 345       array[ipart] = 0.;
 346
 347     Float_t eventtime = 0.;
 348     while(1)
 349       {
 350         eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
 351         if(eventtime > fTimeWindow) break;
 352         array.Set(array.GetSize()+1);
 353         array[array.GetSize()-1] = eventtime;
 354       }
 355
 356     eventtime = 0.;
 357     while(1)
 358       {
 359         eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
 360         if(TMath::Abs(eventtime) > fTimeWindow) break;
 361         array.Set(array.GetSize()+1);
 362         array[array.GetSize()-1] = eventtime;
 363       }
 364
 365     SetNumberParticles(fEventsTime->GetSize());
 366 }
 367
 368 void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
 369                                        Float_t thresh, Float_t etseed, Float_t minet, Float_t r)
 370 {
 371 // Set pycell parameters
 372     fPycellEtaMax    =  etamax;
 373     fPycellNEta      =  neta;
 374     fPycellNPhi      =  nphi;
 375     fPycellThreshold =  thresh;
 376     fPycellEtSeed    =  etseed;
 377     fPycellMinEtJet  =  minet;
 378     fPycellMaxRadius =  r;
 379 }
 380
 381
 382
 383 void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast)
 384 {
 385   // Set a range of event numbers, for which a table
 386   // of generated particle will be printed
 387   fDebugEventFirst = eventFirst;
 388   fDebugEventLast  = eventLast;
 389   if (fDebugEventLast==-1) fDebugEventLast=fDebugEventFirst;
 390 }
 391
 392 void AliGenPythia::Init()
 393 {
 394 // Initialisation
 395
 396     SetMC(AliPythia::Instance());
 397     fPythia=(AliPythia*) fMCEvGen;
 398
 399 //
 400     fParentWeight=1./Float_t(fNpart);
 401 //
 402     if (fWeightPower != 0)
 403       fPythia->SetWeightPower(fWeightPower);
 404     fPythia->SetCKIN(3,fPtHardMin);
 405     fPythia->SetCKIN(4,fPtHardMax);
 406     fPythia->SetCKIN(7,fYHardMin);
 407     fPythia->SetCKIN(8,fYHardMax);
 408
 409     if (fProjectile != "p" || fTarget != "p") fPythia->SetCollisionSystem(fProjectile,fTarget);
 410
 411     if(fUseNuclearPDF)
 412       fPythia->SetNuclei(fAProjectile, fATarget, fNucPdf);
 413     // Fragmentation?
 414     if (fFragmentation) {
 415       fPythia->SetMSTP(111,1);
 416     } else {
 417       fPythia->SetMSTP(111,0);
 418     }
 419
 420
 421 //  initial state radiation
 422     fPythia->SetMSTP(61,fGinit);
 423 //  final state radiation
 424     fPythia->SetMSTP(71,fGfinal);
 425     //color reconnection strength
 426     if(fCRoff==1)fPythia->SetMSTP(95,0);
 427 //  pt - kick
 428     if (fPtKick > 0.) {
 429         fPythia->SetMSTP(91,1);
 430         fPythia->SetPARP(91,fPtKick);
 431         fPythia->SetPARP(93, 4. * fPtKick);
 432     } else {
 433         fPythia->SetMSTP(91,0);
 434     }
 435
 436            if (fReadLHEF) fPythia->OpenFortranFile(97, const_cast<char*>(fkNameLHEF));
 437
 438     if (fReadFromFile) {
 439         fRL  =  AliRunLoader::Open(fkFileName, "Partons");
 440         fRL->LoadKinematics();
 441         fRL->LoadHeader();
 442     } else {
 443         fRL = 0x0;
 444     }
 445  //
 446     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc, fItune);
 447     //  Forward Paramters to the AliPythia object
 448     fDecayer->SetForceDecay(fForceDecay);
 449 // Switch off Heavy Flavors on request
 450     if (fHFoff) {
 451         // Maximum number of quark flavours used in pdf
 452         fPythia->SetMSTP(58, 3);
 453         // Maximum number of flavors that can be used in showers
 454         fPythia->SetMSTJ(45, 3);
 455         // Switch off g->QQbar splitting in decay table
 456         ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
 457     }
 458
 459     fDecayer->Init();
 460
 461
 462 //  Parent and Children Selection
 463     switch (fProcess)
 464     {
 465     case kPyOldUEQ2ordered:
 466     case kPyOldUEQ2ordered2:
 467     case kPyOldPopcorn:
 468       break;
 469     case kPyCharm:
 470     case kPyCharmUnforced:
 471     case kPyCharmPbPbMNR:
 472     case kPyCharmpPbMNR:
 473     case kPyCharmppMNR:
 474     case kPyCharmppMNRwmi:
 475     case kPyCharmPWHG:
 476         fParentSelect[0] =   411;
 477         fParentSelect[1] =   421;
 478         fParentSelect[2] =   431;
 479         fParentSelect[3] =  4122;
 480         fParentSelect[4] =  4232;
 481         fParentSelect[5] =  4132;
 482         fParentSelect[6] =  4332;
 483         fFlavorSelect    =  4;
 484         break;
 485     case kPyD0PbPbMNR:
 486     case kPyD0pPbMNR:
 487     case kPyD0ppMNR:
 488         fParentSelect[0] =   421;
 489         fFlavorSelect    =   4;
 490         break;
 491     case kPyDPlusPbPbMNR:
 492     case kPyDPluspPbMNR:
 493     case kPyDPlusppMNR:
 494         fParentSelect[0] =   411;
 495         fFlavorSelect    =   4;
 496         break;
 497     case kPyDPlusStrangePbPbMNR:
 498     case kPyDPlusStrangepPbMNR:
 499     case kPyDPlusStrangeppMNR:
 500         fParentSelect[0] =   431;
 501         fFlavorSelect    =   4;
 502         break;
 503     case kPyLambdacppMNR:
 504         fParentSelect[0] =  4122;
 505         fFlavorSelect    =   4;
 506         break;
 507     case kPyBeauty:
 508     case kPyBeautyJets:
 509     case kPyBeautyPbPbMNR:
 510     case kPyBeautypPbMNR:
 511     case kPyBeautyppMNR:
 512     case kPyBeautyppMNRwmi:
 513     case kPyBeautyPWHG:
 514         fParentSelect[0]=  511;
 515         fParentSelect[1]=  521;
 516         fParentSelect[2]=  531;
 517         fParentSelect[3]= 5122;
 518         fParentSelect[4]= 5132;
 519         fParentSelect[5]= 5232;
 520         fParentSelect[6]= 5332;
 521         fFlavorSelect   = 5;
 522         break;
 523     case kPyBeautyUnforced:
 524         fParentSelect[0] =  511;
 525         fParentSelect[1] =  521;
 526         fParentSelect[2] =  531;
 527         fParentSelect[3] = 5122;
 528         fParentSelect[4] = 5132;
 529         fParentSelect[5] = 5232;
 530         fParentSelect[6] = 5332;
 531         fFlavorSelect    = 5;
 532         break;
 533     case kPyJpsiChi:
 534     case kPyJpsi:
 535         fParentSelect[0] = 443;
 536         break;
 537     case kPyMbDefault:
 538     case kPyMbAtlasTuneMC09:
 539     case kPyMb:
 540     case kPyMbWithDirectPhoton:
 541     case kPyMbNonDiffr:
 542     case kPyMbMSEL1:
 543     case kPyJets:
 544     case kPyJetsPWHG:
 545     case kPyDirectGamma:
 546     case kPyLhwgMb:
 547         break;
 548     case kPyWPWHG:
 549     case kPyW:
 550     case kPyZ:
 551     case kPyZgamma:
 552     case kPyMBRSingleDiffraction:
 553     case kPyMBRDoubleDiffraction:
 554     case kPyMBRCentralDiffraction:
 555         break;
 556     }
 557 //
 558 //
 559 //  JetFinder for Trigger
 560 //
 561 //  Configure detector (EMCAL like)
 562 //
 563     fPythia->SetPARU(51, fPycellEtaMax);
 564     fPythia->SetMSTU(51, fPycellNEta);
 565     fPythia->SetMSTU(52, fPycellNPhi);
 566 //
 567 //  Configure Jet Finder
 568 //
 569     fPythia->SetPARU(58,  fPycellThreshold);
 570     fPythia->SetPARU(52,  fPycellEtSeed);
 571     fPythia->SetPARU(53,  fPycellMinEtJet);
 572     fPythia->SetPARU(54,  fPycellMaxRadius);
 573     fPythia->SetMSTU(54,  2);
 574 //
 575 //  This counts the total number of calls to Pyevnt() per run.
 576     fTrialsRun = 0;
 577     fQ         = 0.;
 578     fX1        = 0.;
 579     fX2        = 0.;
 580     fNev       = 0 ;
 581 //
 582 //
 583 //
 584     AliGenMC::Init();
 585
 586     // Reset Lorentz boost if demanded
 587     if(!fUseLorentzBoost) {
 588       fDyBoost = 0;
 589       Warning("Init","Demand to discard Lorentz boost.\n");
 590     }
 591 //
 592 //
 593 //
 594     if (fSetNuclei) {
 595       fDyBoost = 0;
 596       Warning("Init","Deprecated function SetNuclei() used (nPDFs + no boost). Use SetProjectile + SetTarget + SetUseNuclearPDF + SetUseLorentzBoost instead.\n");
 597     }
 598     fPythia->SetPARJ(200, 0.0);
 599     fPythia->SetPARJ(199, 0.0);
 600     fPythia->SetPARJ(198, 0.0);
 601     fPythia->SetPARJ(197, 0.0);
 602
 603     if (fQuench == 1) {
 604         fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
 605     }
 606
 607     if(fQuench ==2){fPythia->SetPyquenParameters(fpyquenT,fpyquenTau,fpyquenNf,fpyquenEloss,fpyquenAngle);}
 608
 609     if (fQuench == 3) {
 610         // Nestor's change of the splittings
 611         fPythia->SetPARJ(200, 0.8);
 612         fPythia->SetMSTJ(41, 1);  // QCD radiation only
 613         fPythia->SetMSTJ(42, 2);  // angular ordering
 614         fPythia->SetMSTJ(44, 2);  // option to run alpha_s
 615         fPythia->SetMSTJ(47, 0);  // No correction back to hard scattering element
 616         fPythia->SetMSTJ(50, 0);  // No coherence in first branching
 617         fPythia->SetPARJ(82, 1.); // Cut off for parton showers
 618     } else if (fQuench == 4) {
 619         // Armesto-Cunqueiro-Salgado change of the splittings.
 620         AliFastGlauber* glauber = AliFastGlauber::Instance();
 621         glauber->Init(2);
 622         //read and store transverse almonds corresponding to differnt
 623         //impact parameters.
 624         glauber->SetCentralityClass(0.,0.1);
 625         fPythia->SetPARJ(200, 1.);
 626         fPythia->SetPARJ(198, fQhat);
 627         fPythia->SetPARJ(199, fLength);
 628         fPythia->SetMSTJ(42, 2);  // angular ordering
 629         fPythia->SetMSTJ(44, 2);  // option to run alpha_s
 630         fPythia->SetPARJ(82, 1.); // Cut off for parton showers
 631     }
 632
 633   if ( AliLog::GetDebugLevel("","AliGenPythia") >= 1 ) {
 634     fPythia->Pystat(4);
 635     fPythia->Pystat(5);
 636   }
 637 }
 638
 639 void AliGenPythia::Generate()
 640 {
 641 // Generate one event
 642     if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
 643     fDecayer->ForceDecay();
 644
 645     Double_t polar[3]   =   {0,0,0};
 646     Double_t origin[3]  =   {0,0,0};
 647     Double_t p[4];
 648 //  converts from mm/c to s
 649     const Float_t kconv=0.001/2.999792458e8;
 650 //
 651     Int_t nt=0;
 652     Int_t jev=0;
 653     Int_t j, kf;
 654     fTrials=0;
 655     fEventTime = 0.;
 656
 657
 658
 659     //  Set collision vertex position
 660     if (fVertexSmear == kPerEvent) Vertex();
 661
 662 //  event loop
 663     while(1)
 664     {
 665 //
 666 // Produce event
 667 //
 668 //
 669 // Switch hadronisation off
 670 //
 671         fPythia->SetMSTJ(1, 0);
 672
 673         if (fQuench ==4){
 674             Double_t bimp;
 675             // Quenching comes through medium-modified splitting functions.
 676             AliFastGlauber::Instance()->GetRandomBHard(bimp);
 677             fPythia->SetPARJ(197, bimp);
 678             fImpact = bimp;
 679             fPythia->Qpygin0();
 680         }
 681 //
 682 // Either produce new event or read partons from file
 683 //
 684         if (!fReadFromFile) {
 685             if (!fNewMIS) {
 686                 fPythia->Pyevnt();
 687             } else {
 688                 fPythia->Pyevnw();
 689             }
 690             fNpartons = fPythia->GetN();
 691         } else {
 692             printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
 693             fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
 694             fPythia->SetN(0);
 695             LoadEvent(fRL->Stack(), 0 , 1);
 696             fPythia->Pyedit(21);
 697         }
 698
 699 //
 700 //  Run quenching routine
 701 //
 702         if (fQuench == 1) {
 703             fPythia->Quench();
 704         } else if (fQuench == 2){
 705             fPythia->Pyquen(208., 0, 0.);
 706         } else if (fQuench == 3) {
 707             // Quenching is via multiplicative correction of the splittings
 708         }
 709
 710 //
 711 // Switch hadronisation on
 712 //
 713         if (fHadronisation) {
 714             fPythia->SetMSTJ(1, 1);
 715 //
 716 // .. and perform hadronisation
 717 //      printf("Calling hadronisation %d\n", fPythia->GetN());
 718
 719             if (fPatchOmegaDalitz) {
 720               fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
 721               fPythia->Pyexec();
 722               fPythia->DalitzDecays();
 723               fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
 724             }
 725
 726       else  if (fDecayerExodus) {
 727
 728         fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 0);
 729         fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0);
 730         fPythia->SetMDCY(fPythia->Pycomp(221) ,1, 0);
 731         fPythia->Pyexec();
 732         fPythia->OmegaDalitz();
 733         fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1);
 734         fPythia->PizeroDalitz();
 735         fPythia->PhiDalitz();
 736         fPythia->SetMDCY(fPythia->Pycomp(221) ,1, 1);
 737         fPythia->EtaDalitz();
 738         fPythia->EtaprimeDalitz();
 739         fPythia->SetMDCY(fPythia->Pycomp(22) ,1, 1);
 740         fPythia->RhoDirect();
 741         fPythia->OmegaDirect();
 742         fPythia->PhiDirect();
 743         fPythia->JPsiDirect();
 744       }
 745
 746       fPythia->Pyexec();
 747   }
 748         fTrials++;
 749         fPythia->ImportParticles(&fParticles,"All");
 750
 751         if (TMath::Abs(fDyBoost) > 1.e-4) Boost();
 752         if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle();
 753
 754 //
 755 //
 756 //
 757         Int_t i;
 758
 759         fNprimaries = 0;
 760         Int_t np = fParticles.GetEntriesFast();
 761
 762         if (np == 0) continue;
 763 //
 764
 765 //
 766         Int_t* pParent   = new Int_t[np];
 767         Int_t* pSelected = new Int_t[np];
 768         Int_t* trackIt   = new Int_t[np];
 769         for (i = 0; i < np; i++) {
 770             pParent[i]   = -1;
 771             pSelected[i] =  0;
 772             trackIt[i]   =  0;
 773         }
 774
 775         Int_t nc = 0;        // Total n. of selected particles
 776         Int_t nParents = 0;  // Selected parents
 777         Int_t nTkbles = 0;   // Trackable particles
 778         if (fProcess != kPyMbDefault &&
 779             fProcess != kPyMb &&
 780             fProcess != kPyMbAtlasTuneMC09 &&
 781             fProcess != kPyMbWithDirectPhoton &&
 782             fProcess != kPyJets &&
 783             fProcess != kPyDirectGamma &&
 784             fProcess != kPyMbNonDiffr  &&
 785             fProcess != kPyMbMSEL1     &&
 786             fProcess != kPyW &&
 787             fProcess != kPyZ &&
 788             fProcess != kPyZgamma &&
 789             fProcess != kPyCharmppMNRwmi &&
 790             fProcess != kPyBeautyppMNRwmi &&
 791             fProcess != kPyBeautyJets &&
 792             fProcess != kPyWPWHG &&
 793             fProcess != kPyJetsPWHG &&
 794             fProcess != kPyCharmPWHG &&
 795             fProcess != kPyBeautyPWHG) {
 796
 797             for (i = 0; i < np; i++) {
 798                 TParticle* iparticle = (TParticle *) fParticles.At(i);
 799                 Int_t ks = iparticle->GetStatusCode();
 800                 kf = CheckPDGCode(iparticle->GetPdgCode());
 801 // No initial state partons
 802                 if (ks==21) continue;
 803 //
 804 // Heavy Flavor Selection
 805 //
 806                 // quark ?
 807                 kf = TMath::Abs(kf);
 808                 Int_t kfl = kf;
 809                 // Resonance
 810
 811                 if (kfl > 100000) kfl %= 100000;
 812                 if (kfl > 10000)  kfl %= 10000;
 813                 // meson ?
 814                 if  (kfl > 10) kfl/=100;
 815                 // baryon
 816                 if (kfl > 10) kfl/=10;
 817                 Int_t ipa = iparticle->GetFirstMother()-1;
 818                 Int_t kfMo = 0;
 819 //
 820 // Establish mother daughter relation between heavy quarks and mesons
 821 //
 822                 if (kf >= fFlavorSelect && kf <= 6) {
 823                     Int_t idau = iparticle->GetFirstDaughter() - 1;
 824                     if (idau > -1) {
 825                         TParticle* daughter = (TParticle *) fParticles.At(idau);
 826                         Int_t pdgD = daughter->GetPdgCode();
 827                         if (pdgD == 91 || pdgD == 92) {
 828                             Int_t jmin = daughter->GetFirstDaughter() - 1;
 829                             Int_t jmax = daughter->GetLastDaughter()  - 1;
 830                             for (Int_t jp = jmin; jp <= jmax; jp++)
 831                                 ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
 832                         } // is string or cluster
 833                     } // has daughter
 834                 } // heavy quark
 835
 836
 837                 if (ipa > -1) {
 838                     TParticle *  mother = (TParticle *) fParticles.At(ipa);
 839                     kfMo = TMath::Abs(mother->GetPdgCode());
 840                 }
 841
 842                 // What to keep in Stack?
 843                 Bool_t flavorOK = kFALSE;
 844                 Bool_t selectOK = kFALSE;
 845                 if (fFeedDownOpt) {
 846                     if (kfl >= fFlavorSelect) flavorOK = kTRUE;
 847                 } else {
 848                     if (kfl > fFlavorSelect) {
 849                         nc = -1;
 850                         break;
 851                     }
 852                     if (kfl == fFlavorSelect) flavorOK = kTRUE;
 853                 }
 854                 switch (fStackFillOpt) {
 855                 case kHeavyFlavor:
 856                 case kFlavorSelection:
 857                     selectOK = kTRUE;
 858                     break;
 859                 case kParentSelection:
 860                     if (ParentSelected(kf) || kf <= 10) selectOK = kTRUE;
 861                     break;
 862                 }
 863                 if (flavorOK && selectOK) {
 864 //
 865 // Heavy flavor hadron or quark
 866 //
 867 // Kinematic seletion on final state heavy flavor mesons
 868                     if (ParentSelected(kf) && !KinematicSelection(iparticle, 0))
 869                     {
 870                         continue;
 871                     }
 872                     pSelected[i] = 1;
 873                     if (ParentSelected(kf)) ++nParents; // Update parent count
 874 //                  printf("\n particle (HF)  %d %d %d", i, pSelected[i], kf);
 875                 } else {
 876 // Kinematic seletion on decay products
 877                     if (fCutOnChild && ParentSelected(kfMo) && ChildSelected(kf)
 878                         && !KinematicSelection(iparticle, 1))
 879                     {
 880                         continue;
 881                     }
 882 //
 883 // Decay products
 884 // Select if mother was selected and is not tracked
 885
 886                     if (pSelected[ipa] &&
 887                         !trackIt[ipa]  &&     // mother will be  tracked ?
 888                         kfMo !=  5 &&         // mother is b-quark, don't store fragments
 889                         kfMo !=  4 &&         // mother is c-quark, don't store fragments
 890                         kf   != 92)           // don't store string
 891                     {
 892 //
 893 // Semi-stable or de-selected: diselect decay products:
 894 //
 895 //
 896                         if (pSelected[i] == -1 ||  fDecayer->GetLifetime(kf) > fMaxLifeTime)
 897                         {
 898                             Int_t ipF = iparticle->GetFirstDaughter();
 899                             Int_t ipL = iparticle->GetLastDaughter();
 900                             if (ipF > 0) for (j = ipF-1; j < ipL; j++) pSelected[j] = -1;
 901                         }
 902 //                      printf("\n particle (decay)  %d %d %d", i, pSelected[i], kf);
 903                         pSelected[i] = (pSelected[i] == -1) ? 0 : 1;
 904                     }
 905                 }
 906                 if (pSelected[i] == -1) pSelected[i] = 0;
 907                 if (!pSelected[i]) continue;
 908                 // Count quarks only if you did not include fragmentation
 909                 if (fFragmentation && kf <= 10) continue;
 910
 911                 nc++;
 912 // Decision on tracking
 913                 trackIt[i] = 0;
 914 //
 915 // Track final state particle
 916                 if (ks == 1) trackIt[i] = 1;
 917 // Track semi-stable particles
 918                 if ((ks == 1) || (fDecayer->GetLifetime(kf) > fMaxLifeTime))  trackIt[i] = 1;
 919 // Track particles selected by process if undecayed.
 920                 if (fForceDecay == kNoDecay) {
 921                     if (ParentSelected(kf)) trackIt[i] = 1;
 922                 } else {
 923                     if (ParentSelected(kf)) trackIt[i] = 0;
 924                 }
 925                 if (trackIt[i] == 1) ++nTkbles; // Update trackable counter
 926 //
 927 //
 928
 929             } // particle selection loop
 930             if (nc > 0) {
 931                 for (i = 0; i<np; i++) {
 932                     if (!pSelected[i]) continue;
 933                     TParticle *  iparticle = (TParticle *) fParticles.At(i);
 934                     kf = CheckPDGCode(iparticle->GetPdgCode());
 935                     Int_t ks = iparticle->GetStatusCode();
 936                     p[0] = iparticle->Px();
 937                     p[1] = iparticle->Py();
 938                     p[2] = iparticle->Pz();
 939                     p[3] = iparticle->Energy();
 940
 941                     origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
 942                     origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
 943                     origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
 944
 945                     Float_t tof   = fTime + kconv*iparticle->T();
 946                     Int_t ipa     = iparticle->GetFirstMother()-1;
 947                     Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
 948
 949                     PushTrack(fTrackIt*trackIt[i], iparent, kf,
 950                               p[0], p[1], p[2], p[3],
 951                               origin[0], origin[1], origin[2], tof,
 952                               polar[0], polar[1], polar[2],
 953                               kPPrimary, nt, 1., ks);
 954                     pParent[i] = nt;
 955                     KeepTrack(nt);
 956                     fNprimaries++;
 957                 } //  PushTrack loop
 958             }
 959         } else {
 960             nc = GenerateMB();
 961         } // mb ?
 962
 963         GetSubEventTime();
 964
 965         delete[] pParent;
 966         delete[] pSelected;
 967         delete[] trackIt;
 968
 969         if (nc > 0) {
 970           switch (fCountMode) {
 971           case kCountAll:
 972             // printf(" Count all \n");
 973             jev += nc;
 974             break;
 975           case kCountParents:
 976             // printf(" Count parents \n");
 977             jev += nParents;
 978             break;
 979           case kCountTrackables:
 980             // printf(" Count trackable \n");
 981             jev += nTkbles;
 982             break;
 983           }
 984             if (jev >= fNpart || fNpart == -1) {
 985                 fKineBias=Float_t(fNpart)/Float_t(fTrials);
 986
 987                 fQ  += fPythia->GetVINT(51);
 988                 fX1 += fPythia->GetVINT(41);
 989                 fX2 += fPythia->GetVINT(42);
 990                 fTrialsRun += fTrials;
 991                 fNev++;
 992                 MakeHeader();
 993                 break;
 994             }
 995         }
 996     } // event loop
 997     SetHighWaterMark(nt);
 998 //  adjust weight due to kinematic selection
 999 //    AdjustWeights();
1000 //  get cross-section
1001     fXsection=fPythia->GetPARI(1);
1002 }
1003
1004 Int_t  AliGenPythia::GenerateMB()
1005 {
1006 //
1007 // Min Bias selection and other global selections
1008 //
1009     Int_t i, kf, nt, iparent;
1010     Int_t nc = 0;
1011     Double_t p[4];
1012     Double_t polar[3]   =   {0,0,0};
1013     Double_t origin[3]  =   {0,0,0};
1014 //  converts from mm/c to s
1015     const Float_t kconv=0.001/2.999792458e8;
1016
1017
1018     Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1019
1020
1021
1022     Int_t* pParent = new Int_t[np];
1023     for (i=0; i< np; i++) pParent[i] = -1;
1024     if ((fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG)
1025         && fEtMinJet > 0.) {
1026         TParticle* jet1 = (TParticle *) fParticles.At(6);
1027         TParticle* jet2 = (TParticle *) fParticles.At(7);
1028
1029         if (!jet1 || ! jet2 || !CheckTrigger(jet1, jet2)) {
1030           delete [] pParent;
1031           return 0;
1032         }
1033     }
1034
1035
1036   // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
1037   // implemented primaryly for kPyJets, but extended to any kind of process.
1038   if ((fFragPhotonInCalo || fPi0InCalo || fEtaInCalo || fEleInCalo || fHadronInCalo || fDecayPhotonInCalo) &&
1039       (fCheckPHOS || fCheckEMCAL || fCheckBarrel) ) {
1040     Bool_t ok = TriggerOnSelectedParticles(np);
1041
1042     if(!ok) {
1043       delete[] pParent;
1044       return 0;
1045     }
1046   }
1047
1048     // Check for diffraction
1049     if(fkTuneForDiff) {
1050       if(fItune==320 && ( (TMath::Abs(fEnergyCMS - 900) < 1) || (TMath::Abs(fEnergyCMS - 2760) < 1) || (TMath::Abs(fEnergyCMS - 7000) < 1)) ) {
1051         if(!CheckDiffraction()) {
1052           delete [] pParent;
1053           return 0;
1054         }
1055       }
1056     }
1057
1058     // Check for minimum multiplicity
1059     if (fTriggerMultiplicity > 0) {
1060       Int_t multiplicity = 0;
1061       for (i = 0; i < np; i++) {
1062         TParticle *  iparticle = (TParticle *) fParticles.At(i);
1063
1064         Int_t statusCode = iparticle->GetStatusCode();
1065
1066         // Initial state particle
1067         if (statusCode != 1)
1068           continue;
1069         // eta cut
1070         if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
1071           continue;
1072         //multiplicity check for a given eta range
1073         if ((fTriggerMultiplicityEtaMin != fTriggerMultiplicityEtaMax) &&
1074             (iparticle->Eta() < fTriggerMultiplicityEtaMin || iparticle->Eta() > fTriggerMultiplicityEtaMax))
1075           continue;
1076         // pt cut
1077         if (iparticle->Pt() < fTriggerMultiplicityPtMin)
1078             continue;
1079
1080         TParticlePDG* pdgPart = iparticle->GetPDG();
1081         if (pdgPart && pdgPart->Charge() == 0)
1082           continue;
1083
1084         ++multiplicity;
1085       }
1086
1087       if (multiplicity < fTriggerMultiplicity) {
1088         delete [] pParent;
1089         return 0;
1090       }
1091       Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity);
1092     }
1093
1094
1095     if (fTriggerParticle) {
1096         Bool_t triggered = kFALSE;
1097         for (i = 0; i < np; i++) {
1098             TParticle *  iparticle = (TParticle *) fParticles.At(i);
1099             kf = CheckPDGCode(iparticle->GetPdgCode());
1100             if (kf != fTriggerParticle) continue;
1101             if (iparticle->Pt() == 0.) continue;
1102             if (TMath::Abs(iparticle->Y()) > fTriggerY) continue;
1103             if (fTriggerEtaMin == fTriggerEta) {
1104               if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
1105             } else {
1106               if (iparticle->Eta() > fTriggerEta || iparticle->Eta() < fTriggerEtaMin) continue;
1107             }
1108             if ( iparticle->Pt() > fTriggerMaxPt || iparticle->Pt() < fTriggerMinPt ) continue;
1109             triggered = kTRUE;
1110             break;
1111         }
1112         if (!triggered) {
1113           delete [] pParent;
1114           return 0;
1115         }
1116     }
1117
1118
1119     // Check if there is a ccbar or bbbar pair with at least one of the two
1120     // in fYMin < y < fYMax
1121
1122     if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) {
1123       TParticle *partCheck;
1124       TParticle *mother;
1125       Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
1126       Bool_t  theChild=kFALSE;
1127       Bool_t  triggered=kFALSE;
1128       Float_t y;
1129       Int_t   pdg, mpdg, mpdgUpperFamily;
1130       for(i = 0; i < np; i++) {
1131         partCheck = (TParticle*)fParticles.At(i);
1132         pdg = partCheck->GetPdgCode();
1133         if(TMath::Abs(pdg) == fFlavorSelect) { // quark
1134           if(pdg>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
1135           y = 0.5*TMath::Log((partCheck->Energy()+partCheck->Pz()+1.e-13)/
1136                              (partCheck->Energy()-partCheck->Pz()+1.e-13));
1137           if(fUseYCutHQ && y>fYMinHQ && y<fYMaxHQ) inYcut=kTRUE;
1138           if(!fUseYCutHQ && y>fYMin && y<fYMax) inYcut=kTRUE;
1139         }
1140         if(fTriggerParticle) {
1141           if(TMath::Abs(pdg)==fTriggerParticle) triggered=kTRUE;
1142         }
1143         if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
1144           Int_t mi = partCheck->GetFirstMother() - 1;
1145           if(mi<0) continue;
1146           mother = (TParticle*)fParticles.At(mi);
1147           mpdg=TMath::Abs(mother->GetPdgCode());
1148           mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
1149           if ( ParentSelected(mpdg) ||
1150               (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
1151             if (KinematicSelection(partCheck,1)) {
1152               theChild=kTRUE;
1153             }
1154           }
1155         }
1156       }
1157       if (!theQ || !theQbar || !inYcut) { // one of the c/b conditions not satisfied
1158         delete[] pParent;
1159         return 0;
1160       }
1161       if (fCutOnChild && !theChild) { // one of the child conditions not satisfied
1162         delete[] pParent;
1163         return 0;
1164       }
1165       if(fTriggerParticle && !triggered) { // particle requested is not produced
1166         delete[] pParent;
1167         return 0;
1168       }
1169
1170     }
1171
1172     //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
1173     if ( (
1174           fProcess == kPyWPWHG ||
1175           fProcess == kPyW ||
1176           fProcess == kPyZ ||
1177           fProcess == kPyZgamma ||
1178           fProcess == kPyMbDefault ||
1179           fProcess == kPyMb ||
1180           fProcess == kPyMbAtlasTuneMC09 ||
1181           fProcess == kPyMbWithDirectPhoton ||
1182           fProcess == kPyMbNonDiffr)
1183          && (fCutOnChild == 1) ) {
1184       if ( !CheckKinematicsOnChild() ) {
1185         delete[] pParent;
1186         return 0;
1187       }
1188     }
1189
1190
1191     for (i = 0; i < np; i++) {
1192         Int_t trackIt = 0;
1193         TParticle *  iparticle = (TParticle *) fParticles.At(i);
1194         kf = CheckPDGCode(iparticle->GetPdgCode());
1195         Int_t ks = iparticle->GetStatusCode();
1196         Int_t km = iparticle->GetFirstMother();
1197         if (
1198             (((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) || (ks !=1)) && ((fStackFillOpt != kHeavyFlavor) || IsFromHeavyFlavor(i))) ||
1199             ((fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG) && ks == 21 && km == 0 && i>1)
1200             )
1201           {
1202             nc++;
1203             if (ks == 1) trackIt = 1;
1204             Int_t ipa = iparticle->GetFirstMother()-1;
1205
1206             iparent = (ipa > -1) ? pParent[ipa] : -1;
1207
1208 //
1209 // store track information
1210             p[0] = iparticle->Px();
1211             p[1] = iparticle->Py();
1212             p[2] = iparticle->Pz();
1213             p[3] = iparticle->Energy();
1214
1215
1216             origin[0] = fVertex[0]+iparticle->Vx()/10; // [cm]
1217             origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
1218             origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
1219
1220             Float_t tof = fTime + fEventTime + kconv * iparticle->T();
1221
1222             PushTrack(fTrackIt*trackIt, iparent, kf,
1223                       p[0], p[1], p[2], p[3],
1224                       origin[0], origin[1], origin[2], tof,
1225                       polar[0], polar[1], polar[2],
1226                       kPPrimary, nt, 1., ks);
1227             fNprimaries++;
1228             KeepTrack(nt);
1229             pParent[i] = nt;
1230             SetHighWaterMark(nt);
1231
1232         } // select particle
1233     } // particle loop
1234
1235     delete[] pParent;
1236
1237     return 1;
1238 }
1239
1240
1241 void AliGenPythia::FinishRun()
1242 {
1243 // Print x-section summary
1244     fPythia->Pystat(1);
1245
1246     if (fNev > 0.) {
1247         fQ  /= fNev;
1248         fX1 /= fNev;
1249         fX2 /= fNev;
1250     }
1251
1252     printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
1253     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
1254 }
1255
1256 void AliGenPythia::AdjustWeights() const
1257 {
1258 // Adjust the weights after generation of all events
1259 //
1260     if (gAlice) {
1261         TParticle *part;
1262         Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
1263         for (Int_t i=0; i<ntrack; i++) {
1264             part= gAlice->GetMCApp()->Particle(i);
1265             part->SetWeight(part->GetWeight()*fKineBias);
1266         }
1267     }
1268 }
1269
1270 void AliGenPythia::SetNuclei(Int_t a1, Int_t a2, Int_t pdfset)
1271 {
1272 // Treat protons as inside nuclei with mass numbers a1 and a2
1273
1274     fAProjectile   = a1;
1275     fATarget       = a2;
1276     fNucPdf        = pdfset;  // 0 EKS98 9 EPS09LO 19 EPS09NLO
1277     fUseNuclearPDF = kTRUE;
1278     fSetNuclei     = kTRUE;
1279 }
1280
1281
1282 void AliGenPythia::MakeHeader()
1283 {
1284 //
1285 // Make header for the simulated event
1286 //
1287   if (gAlice) {
1288     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
1289         gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
1290   }
1291
1292 // Builds the event header, to be called after each event
1293     if (fHeader) delete fHeader;
1294     fHeader = new AliGenPythiaEventHeader("Pythia");
1295     fHeader->SetTitle(GetTitle());
1296 //
1297 // Event type
1298     if(fProcDiff>0){
1299       //      if(fProcDiff == 92 || fProcDiff == 93) printf("\n\n\n\n\n");
1300       //      printf("fPythia->GetMSTI(1) = %d   fProcDiff = %d\n",fPythia->GetMSTI(1), fProcDiff);
1301     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fProcDiff);
1302     }
1303     else
1304     ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
1305 //
1306 // Number of trials
1307     ((AliGenPythiaEventHeader*) fHeader)->SetTrials(fTrials);
1308 //
1309 // Event Vertex
1310     fHeader->SetPrimaryVertex(fVertex);
1311     fHeader->SetInteractionTime(fTime+fEventTime);
1312 //
1313 // Number of primaries
1314     fHeader->SetNProduced(fNprimaries);
1315 //
1316 // Event weight
1317     fHeader->SetEventWeight(fPythia->GetVINT(97));
1318 //
1319 // Jets that have triggered
1320
1321     //Need to store jets for b-jet studies too!
1322     if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG)
1323     {
1324         Int_t ntrig, njet;
1325         Float_t jets[4][10];
1326         GetJets(njet, ntrig, jets);
1327
1328
1329         for (Int_t i = 0; i < ntrig; i++) {
1330             ((AliGenPythiaEventHeader*) fHeader)->AddJet(jets[0][i], jets[1][i], jets[2][i],
1331                                                         jets[3][i]);
1332         }
1333     }
1334 //
1335 // Copy relevant information from external header, if present.
1336 //
1337     Float_t uqJet[4];
1338
1339     if (fRL) {
1340         AliGenPythiaEventHeader* exHeader = (AliGenPythiaEventHeader*) (fRL->GetHeader()->GenEventHeader());
1341         for (Int_t i = 0; i < exHeader->NTriggerJets(); i++)
1342         {
1343             printf("Adding Jet %d %d \n", i,  exHeader->NTriggerJets());
1344
1345
1346             exHeader->TriggerJet(i, uqJet);
1347             ((AliGenPythiaEventHeader*) fHeader)->AddUQJet(uqJet[0], uqJet[1], uqJet[2], uqJet[3]);
1348         }
1349     }
1350 //
1351 // Store quenching parameters
1352 //
1353     if (fQuench){
1354         Double_t z[4] = {0.};
1355         Double_t xp = 0.;
1356         Double_t yp = 0.;
1357
1358         if (fQuench == 1) {
1359             // Pythia::Quench()
1360             fPythia->GetQuenchingParameters(xp, yp, z);
1361         } else if (fQuench == 2){
1362             // Pyquen
1363             Double_t r1 = PARIMP.rb1;
1364             Double_t r2 = PARIMP.rb2;
1365             Double_t b  = PARIMP.b1;
1366             Double_t r   = 0.5 * TMath::Sqrt(2. * (r1 * r1 + r2 * r2) - b * b);
1367             Double_t phi = PARIMP.psib1;
1368             xp = r * TMath::Cos(phi);
1369             yp = r * TMath::Sin(phi);
1370
1371         } else if (fQuench == 4) {
1372             // QPythia
1373             Double_t xy[2];
1374             Double_t i0i1[2];
1375             AliFastGlauber::Instance()->GetSavedXY(xy);
1376             AliFastGlauber::Instance()->GetSavedI0I1(i0i1);
1377             xp = xy[0];
1378             yp = xy[1];
1379             ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact);
1380         }
1381
1382             ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
1383             ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
1384     }
1385 //
1386 // Store pt^hard and cross-section
1387     ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
1388     ((AliGenPythiaEventHeader*) fHeader)->SetXsection(fPythia->GetPARI(1));
1389
1390 //
1391 // Store Event Weight
1392     ((AliGenPythiaEventHeader*) fHeader)->SetEventWeight(fPythia->GetPARI(7)*fPythia->GetPARI(10));
1393     // PARI(7) is 1 or -1, for weighted generation with accept/reject, e.g. POWHEG
1394     // PARI(10) is a weight associated with reweighted generation, using Pyevwt
1395 //
1396 //  Pass header
1397 //
1398     AddHeader(fHeader);
1399     fHeader = 0x0;
1400 }
1401
1402 Bool_t AliGenPythia::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
1403 {
1404 // Check the kinematic trigger condition
1405 //
1406     Double_t eta[2];
1407     eta[0] = jet1->Eta();
1408     eta[1] = jet2->Eta();
1409     Double_t phi[2];
1410     phi[0] = jet1->Phi();
1411     phi[1] = jet2->Phi();
1412     Int_t    pdg[2];
1413     pdg[0] = jet1->GetPdgCode();
1414     pdg[1] = jet2->GetPdgCode();
1415     Bool_t   triggered = kFALSE;
1416
1417     if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess ==  kPyBeautyppMNRwmi || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG) {
1418         Int_t njets = 0;
1419         Int_t ntrig = 0;
1420         Float_t jets[4][10];
1421 //
1422 // Use Pythia clustering on parton level to determine jet axis
1423 //
1424         GetJets(njets, ntrig, jets);
1425
1426         if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
1427 //
1428     } else {
1429         Int_t ij = 0;
1430         Int_t ig = 1;
1431         if (pdg[0] == kGamma) {
1432             ij = 1;
1433             ig = 0;
1434         }
1435         //Check eta range first...
1436         if ((eta[ij] < fEtaMaxJet   && eta[ij] > fEtaMinJet) &&
1437             (eta[ig] < fEtaMaxGamma && eta[ig] > fEtaMinGamma))
1438         {
1439             //Eta is okay, now check phi range
1440             if ((phi[ij] < fPhiMaxJet   && phi[ij] > fPhiMinJet) &&
1441                 (phi[ig] < fPhiMaxGamma && phi[ig] > fPhiMinGamma))
1442             {
1443                 triggered = kTRUE;
1444             }
1445         }
1446     }
1447     return triggered;
1448 }
1449
1450
1451
1452 Bool_t AliGenPythia::CheckKinematicsOnChild(){
1453 //
1454 //Checking Kinematics on Child (status code 1, particle code ?, kin cuts
1455 //
1456     Bool_t checking = kFALSE;
1457     Int_t j, kcode, ks, km;
1458     Int_t nPartAcc = 0; //number of particles in the acceptance range
1459     Int_t numberOfAcceptedParticles = 1;
1460     if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
1461     Int_t npart = fParticles.GetEntriesFast();
1462
1463     for (j = 0; j<npart; j++) {
1464         TParticle *  jparticle = (TParticle *) fParticles.At(j);
1465         kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
1466         ks = jparticle->GetStatusCode();
1467         km = jparticle->GetFirstMother();
1468
1469         if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
1470             nPartAcc++;
1471         }
1472         if( numberOfAcceptedParticles <= nPartAcc){
1473           checking = kTRUE;
1474           break;
1475         }
1476     }
1477
1478     return checking;
1479 }
1480
1481 void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
1482 {
1483   //
1484   // Load event into Pythia Common Block
1485   //
1486
1487   Int_t npart = stack -> GetNprimary();
1488   Int_t n0 = 0;
1489
1490   if (!flag) {
1491     (fPythia->GetPyjets())->N = npart;
1492   } else {
1493     n0 = (fPythia->GetPyjets())->N;
1494     (fPythia->GetPyjets())->N = n0 + npart;
1495   }
1496
1497
1498   for (Int_t part = 0; part < npart; part++) {
1499     TParticle *mPart = stack->Particle(part);
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 void  AliGenPythia::LoadEvent(const TObjArray* stack, Int_t flag, Int_t reHadr)
1536 {
1537   //
1538   // Load event into Pythia Common Block
1539   //
1540
1541   Int_t npart = stack -> GetEntries();
1542   Int_t n0 = 0;
1543
1544   if (!flag) {
1545     (fPythia->GetPyjets())->N = npart;
1546   } else {
1547     n0 = (fPythia->GetPyjets())->N;
1548     (fPythia->GetPyjets())->N = n0 + npart;
1549   }
1550
1551
1552   for (Int_t part = 0; part < npart; part++) {
1553     TParticle *mPart = dynamic_cast<TParticle *>(stack->At(part));
1554     if (!mPart) continue;
1555
1556     Int_t kf     =  mPart->GetPdgCode();
1557     Int_t ks     =  mPart->GetStatusCode();
1558     Int_t idf    =  mPart->GetFirstDaughter();
1559     Int_t idl    =  mPart->GetLastDaughter();
1560
1561     if (reHadr) {
1562         if (ks == 11 || ks == 12) {
1563             ks  -= 10;
1564             idf  = -1;
1565             idl  = -1;
1566         }
1567     }
1568
1569     Float_t px = mPart->Px();
1570     Float_t py = mPart->Py();
1571     Float_t pz = mPart->Pz();
1572     Float_t e  = mPart->Energy();
1573     Float_t m  = mPart->GetCalcMass();
1574
1575
1576     (fPythia->GetPyjets())->P[0][part+n0] = px;
1577     (fPythia->GetPyjets())->P[1][part+n0] = py;
1578     (fPythia->GetPyjets())->P[2][part+n0] = pz;
1579     (fPythia->GetPyjets())->P[3][part+n0] = e;
1580     (fPythia->GetPyjets())->P[4][part+n0] = m;
1581
1582     (fPythia->GetPyjets())->K[1][part+n0] = kf;
1583     (fPythia->GetPyjets())->K[0][part+n0] = ks;
1584     (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
1585     (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
1586     (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
1587   }
1588 }
1589
1590
1591 void AliGenPythia::RecJetsUA1(Int_t& njets, Float_t jets [4][50])
1592 {
1593 //
1594 //  Calls the Pythia jet finding algorithm to find jets in the current event
1595 //
1596 //
1597 //
1598 //  Save jets
1599     Int_t n     = fPythia->GetN();
1600
1601 //
1602 //  Run Jet Finder
1603     fPythia->Pycell(njets);
1604     Int_t i;
1605     for (i = 0; i < njets; i++) {
1606         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1607         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1608         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1609         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1610
1611         jets[0][i] = px;
1612         jets[1][i] = py;
1613         jets[2][i] = pz;
1614         jets[3][i] = e;
1615     }
1616 }
1617
1618
1619
1620 void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
1621 {
1622 //
1623 //  Calls the Pythia clustering algorithm to find jets in the current event
1624 //
1625     Int_t n     = fPythia->GetN();
1626     nJets       = 0;
1627     nJetsTrig   = 0;
1628     if (fJetReconstruction == kCluster) {
1629 //
1630 //  Configure cluster algorithm
1631 //
1632         fPythia->SetPARU(43, 2.);
1633         fPythia->SetMSTU(41, 1);
1634 //
1635 //  Call cluster algorithm
1636 //
1637         fPythia->Pyclus(nJets);
1638 //
1639 //  Loading jets from common block
1640 //
1641     } else {
1642
1643 //
1644 //  Run Jet Finder
1645         fPythia->Pycell(nJets);
1646     }
1647
1648     Int_t i;
1649     for (i = 0; i < nJets; i++) {
1650         Float_t px    = (fPythia->GetPyjets())->P[0][n+i];
1651         Float_t py    = (fPythia->GetPyjets())->P[1][n+i];
1652         Float_t pz    = (fPythia->GetPyjets())->P[2][n+i];
1653         Float_t e     = (fPythia->GetPyjets())->P[3][n+i];
1654         Float_t pt    = TMath::Sqrt(px * px + py * py);
1655         Float_t phi   = TMath::Pi() + TMath::ATan2(-py, -px);
1656         Float_t theta = TMath::ATan2(pt,pz);
1657         Float_t et    = e * TMath::Sin(theta);
1658         Float_t eta   = -TMath::Log(TMath::Tan(theta / 2.));
1659         if (
1660             eta > fEtaMinJet && eta < fEtaMaxJet &&
1661             phi > fPhiMinJet && phi < fPhiMaxJet &&
1662             et  > fEtMinJet  && et  < fEtMaxJet
1663             )
1664         {
1665             jets[0][nJetsTrig] = px;
1666             jets[1][nJetsTrig] = py;
1667             jets[2][nJetsTrig] = pz;
1668             jets[3][nJetsTrig] = e;
1669             nJetsTrig++;
1670 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1671         } else {
1672 //          printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg);
1673         }
1674     }
1675 }
1676
1677 void AliGenPythia::GetSubEventTime()
1678 {
1679   // Calculates time of the next subevent
1680   fEventTime = 0.;
1681   if (fEventsTime) {
1682     TArrayF &array = *fEventsTime;
1683     fEventTime = array[fCurSubEvent++];
1684   }
1685   //  printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
1686   return;
1687 }
1688
1689 Bool_t AliGenPythia::IsInBarrel(Float_t eta) const
1690 {
1691   // Is particle in Central Barrel acceptance?
1692   // etamin=-etamax
1693   if( eta < fTriggerEta  )
1694     return kTRUE;
1695   else
1696     return kFALSE;
1697 }
1698
1699 Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta) const
1700 {
1701   // Is particle in EMCAL acceptance?
1702   // phi in degrees, etamin=-etamax
1703   if(phi > fEMCALMinPhi  && phi < fEMCALMaxPhi &&
1704      eta < fEMCALEta  )
1705     return kTRUE;
1706   else
1707     return kFALSE;
1708 }
1709
1710 Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta, Int_t iparticle)
1711 {
1712   // Is particle in PHOS acceptance?
1713   // Acceptance slightly larger considered.
1714   // phi in degrees, etamin=-etamax
1715   // iparticle is the index of the particle to be checked, for PHOS rotation case
1716
1717   if(phi > fPHOSMinPhi  && phi < fPHOSMaxPhi &&
1718      eta < fPHOSEta  )
1719     return kTRUE;
1720   else
1721   {
1722     if( fCheckPHOSeta && eta < fPHOSEta) fPHOSRotateCandidate = iparticle;
1723
1724     return kFALSE;
1725   }
1726 }
1727
1728 void AliGenPythia::RotatePhi(Bool_t& okdd)
1729 {
1730   //Rotate event in phi to enhance events in PHOS acceptance
1731
1732   if(fPHOSRotateCandidate < 0) return ;
1733
1734   //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
1735   Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
1736   Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
1737   Double_t phiPHOS = (AliPythiaRndm::GetPythiaRandom())->Uniform(phiPHOSmin,phiPHOSmax);
1738
1739   //calculate deltaphi
1740   TParticle* ph = (TParticle *) fParticles.At(fPHOSRotateCandidate);
1741   Double_t phphi = ph->Phi();
1742   Double_t deltaphi = phiPHOS - phphi;
1743
1744
1745
1746   //loop for all particles and produce the phi rotation
1747   Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1748   Double_t oldphi, newphi;
1749   Double_t newVx, newVy, r, vZ, time;
1750   Double_t newPx, newPy, pt, pz, e;
1751   for(Int_t i=0; i< np; i++) {
1752     TParticle* iparticle = (TParticle *) fParticles.At(i);
1753     oldphi = iparticle->Phi();
1754     newphi = oldphi + deltaphi;
1755     if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
1756     if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
1757
1758     r = iparticle->R();
1759     newVx = r * TMath::Cos(newphi);
1760     newVy = r * TMath::Sin(newphi);
1761     vZ   = iparticle->Vz(); // don't transform
1762     time = iparticle->T(); // don't transform
1763
1764     pt = iparticle->Pt();
1765     newPx = pt * TMath::Cos(newphi);
1766     newPy = pt * TMath::Sin(newphi);
1767     pz = iparticle->Pz(); // don't transform
1768     e  = iparticle->Energy(); // don't transform
1769
1770     // apply rotation
1771     iparticle->SetProductionVertex(newVx, newVy, vZ, time);
1772     iparticle->SetMomentum(newPx, newPy, pz, e);
1773
1774   } //end particle loop
1775
1776   // now let's check that we put correctly the candidate photon in PHOS
1777   Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
1778   Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
1779   if(IsInPHOS(phi,eta,-1))
1780     okdd = kTRUE;
1781
1782   // reset the value for next event
1783   fPHOSRotateCandidate = -1;
1784
1785 }
1786
1787
1788 Bool_t AliGenPythia::CheckDiffraction()
1789 {
1790   // use this method only with Perugia-0 tune!
1791
1792   //  printf("AAA\n");
1793
1794    Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
1795
1796    Int_t iPart1=-1;
1797    Int_t iPart2=-1;
1798
1799    Double_t y1 = 1e10;
1800    Double_t y2 = -1e10;
1801
1802   const Int_t kNstable=20;
1803   const Int_t pdgStable[20] = {
1804     22,             // Photon
1805     11,             // Electron
1806     12,             // Electron Neutrino
1807     13,             // Muon
1808     14,             // Muon Neutrino
1809     15,             // Tau
1810     16,             // Tau Neutrino
1811     211,            // Pion
1812     321,            // Kaon
1813     311,            // K0
1814     130,            // K0s
1815     310,            // K0l
1816     2212,           // Proton
1817     2112,           // Neutron
1818     3122,           // Lambda_0
1819     3112,           // Sigma Minus
1820     3222,           // Sigma Plus
1821     3312,           // Xsi Minus
1822     3322,           // Xsi0
1823     3334            // Omega
1824   };
1825
1826      for (Int_t i = 0; i < np; i++) {
1827         TParticle *  part = (TParticle *) fParticles.At(i);
1828
1829         Int_t statusCode = part->GetStatusCode();
1830
1831         // Initial state particle
1832         if (statusCode != 1)
1833           continue;
1834
1835         Int_t pdg = TMath::Abs(part->GetPdgCode());
1836         Bool_t isStable = kFALSE;
1837         for (Int_t i1 = 0; i1 < kNstable; i1++) {
1838           if (pdg == pdgStable[i1]) {
1839             isStable = kTRUE;
1840             break;
1841           }
1842         }
1843         if(!isStable)
1844           continue;
1845
1846         Double_t y = part->Y();
1847
1848         if (y < y1)
1849           {
1850             y1 = y;
1851             iPart1 = i;
1852           }
1853         if (y > y2)
1854         {
1855           y2 = y;
1856           iPart2 = i;
1857         }
1858      }
1859
1860      if(iPart1<0 || iPart2<0) return kFALSE;
1861
1862      y1=TMath::Abs(y1);
1863      y2=TMath::Abs(y2);
1864
1865      TParticle *  part1 = (TParticle *) fParticles.At(iPart1);
1866      TParticle *  part2 = (TParticle *) fParticles.At(iPart2);
1867
1868      Int_t pdg1 = part1->GetPdgCode();
1869      Int_t pdg2 = part2->GetPdgCode();
1870
1871
1872      Int_t iPart = -1;
1873      if (pdg1 == 2212 && pdg2 == 2212)
1874        {
1875          if(y1 > y2)
1876            iPart = iPart1;
1877          else if(y1 < y2)
1878            iPart = iPart2;
1879          else {
1880            iPart = iPart1;
1881            if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)>0.5) iPart = iPart2;
1882          }
1883        }
1884      else if (pdg1 == 2212)
1885        iPart = iPart1;
1886      else if (pdg2 == 2212)
1887        iPart = iPart2;
1888
1889
1890
1891
1892
1893      Double_t M=-1.;
1894      if(iPart>0) {
1895        TParticle *  part = (TParticle *) fParticles.At(iPart);
1896        Double_t E= part->Energy();
1897        Double_t P= part->P();
1898        Double_t M2 = (fEnergyCMS-E-P)*(fEnergyCMS-E+P);
1899        if(M2<0)  return kFALSE;
1900        M= TMath::Sqrt(M2);
1901      }
1902
1903      Double_t Mmin, Mmax, wSD, wDD, wND;
1904      if(!GetWeightsDiffraction(M, Mmin, Mmax, wSD, wDD, wND)) return kFALSE;
1905
1906      if(M>-1 && M<Mmin) return kFALSE;
1907      if(M>Mmax) M=-1;
1908
1909      Int_t procType=fPythia->GetMSTI(1);
1910      Int_t proc0=2;
1911      if(procType== 94) proc0=1;
1912      if(procType== 92 || procType== 93) proc0=0;
1913
1914      Int_t proc=2;
1915      if(M>0) proc=0;
1916      else if(proc0==1) proc=1;
1917
1918      if(proc==1 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wDD) return kFALSE;
1919      if(proc==2 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wND) return kFALSE;
1920
1921
1922     //     if(proc==1 || proc==2) return kFALSE;
1923
1924      if(proc!=0) {
1925        if(proc0!=0) fProcDiff = procType;
1926        else       fProcDiff = 95;
1927        return kTRUE;
1928      }
1929
1930     if(wSD<0)  AliError("wSD<0 ! \n");
1931
1932     if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)> wSD) return kFALSE;
1933
1934     //    printf("iPart = %d\n", iPart);
1935
1936     if(iPart==iPart1) fProcDiff=93;
1937     else if(iPart==iPart2) fProcDiff=92;
1938     else {
1939       printf("EROOR:  iPart!=iPart1 && iPart!=iPart2\n");
1940
1941     }
1942
1943     return kTRUE;
1944 }
1945
1946
1947
1948 Bool_t AliGenPythia::GetWeightsDiffraction(Double_t M, Double_t &Mmin, Double_t &Mmax,
1949                                                        Double_t &wSD, Double_t &wDD, Double_t &wND)
1950 {
1951
1952   // 900 GeV
1953   if(TMath::Abs(fEnergyCMS-900)<1 ){
1954
1955 const Int_t nbin=400;
1956 Double_t bin[]={
1957 1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
1958 4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
1959 7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
1960 10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
1961 13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
1962 15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
1963 18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
1964 21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
1965 24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
1966 27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
1967 30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
1968 33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
1969 36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
1970 39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
1971 42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
1972 45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
1973 48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
1974 51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
1975 54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
1976 57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
1977 60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
1978 63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
1979 66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
1980 69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
1981 72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
1982 75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
1983 78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
1984 81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
1985 84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
1986 87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
1987 90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
1988 93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
1989 96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
1990 99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
1991 102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
1992 105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
1993 108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
1994 111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
1995 114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
1996 117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
1997 120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
1998 123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
1999 126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2000 129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2001 132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2002 135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2003 138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2004 141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2005 144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2006 147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2007 150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2008 153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2009 156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2010 159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2011 162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2012 165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2013 168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2014 171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2015 174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2016 177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2017 180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2018 183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2019 186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2020 189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2021 192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2022 195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2023 198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2024 Double_t w[]={
2025 1.000000, 0.643457, 0.645609, 0.648347, 0.604563, 0.605002,
2026 0.602819, 0.611473, 0.576412, 0.562354, 0.550216, 0.529285,
2027 0.534558, 0.534364, 0.530358, 0.518475, 0.489253, 0.469754,
2028 0.469825, 0.450513, 0.455849, 0.435312, 0.437210, 0.456686,
2029 0.413577, 0.427093, 0.426894, 0.418834, 0.409475, 0.388483,
2030 0.412208, 0.388912, 0.389611, 0.382572, 0.389220, 0.370964,
2031 0.380463, 0.370873, 0.363701, 0.369363, 0.357361, 0.365759,
2032 0.348566, 0.337062, 0.348190, 0.332330, 0.359001, 0.335836,
2033 0.339154, 0.335599, 0.336035, 0.335204, 0.353440, 0.337836,
2034 0.333874, 0.307120, 0.294963, 0.324978, 0.313359, 0.317688,
2035 0.323758, 0.319304, 0.335317, 0.301765, 0.317257, 0.356331,
2036 0.323119, 0.297732, 0.303188, 0.326102, 0.316467, 0.294728,
2037 0.308135, 0.288857, 0.325692, 0.312493, 0.291100, 0.325921,
2038 0.313317, 0.295980, 0.308481, 0.328380, 0.313081, 0.296763,
2039 0.295431, 0.317325, 0.320462, 0.286918, 0.316035, 0.335208,
2040 0.283193, 0.333945, 0.292534, 0.294164, 0.330680, 0.296992,
2041 0.285509, 0.317260, 0.298269, 0.311299, 0.312129, 0.286822,
2042 0.287442, 0.319139, 0.283314, 0.318454, 0.297727, 0.301597,
2043 0.282483, 0.294792, 0.305569, 0.290957, 0.297817, 0.282908,
2044 0.272401, 0.305584, 0.300220, 0.297020, 0.298781, 0.278008,
2045 0.277727, 0.323777, 0.287419, 0.342074, 0.287259, 0.303658,
2046 0.302668, 0.279622, 0.280586, 0.313630, 0.276068, 0.257051,
2047 0.309996, 0.265534, 0.297138, 0.281738, 0.294610, 0.292882,
2048 0.286860, 0.312686, 0.293244, 0.293744, 0.271375, 0.278734,
2049 0.280308, 0.304739, 0.287907, 0.285261, 0.311180, 0.313476,
2050 0.289660, 0.289282, 0.319505, 0.271285, 0.272008, 0.289245,
2051 0.281038, 0.285284, 0.295836, 0.281416, 0.283501, 0.295417,
2052 0.304372, 0.297764, 0.291378, 0.321530, 0.315604, 0.329507,
2053 0.282609, 0.275576, 0.283721, 0.311714, 0.283760, 0.273188,
2054 0.312193, 0.264347, 0.281532, 0.301226, 0.281718, 0.336408,
2055 0.283157, 0.332010, 0.289974, 0.290256, 0.301569, 0.332228,
2056 0.288282, 0.326339, 0.313653, 0.300361, 0.289470, 0.264830,
2057 0.298659, 0.272359, 0.278878, 0.306001, 0.328168, 0.294991,
2058 0.327737, 0.278056, 0.302435, 0.284183, 0.279270, 0.307279,
2059 0.307917, 0.315196, 0.283803, 0.313333, 0.315730, 0.304818,
2060 0.307171, 0.295223, 0.333741, 0.346911, 0.310143, 0.336686,
2061 0.275459, 0.334781, 0.295405, 0.275816, 0.301678, 0.327242,
2062 0.320717, 0.309230, 0.292145, 0.294489, 0.305088, 0.300969,
2063 0.277438, 0.326159, 0.297065, 0.301177, 0.303843, 0.275382,
2064 0.304019, 0.284166, 0.289610, 0.331611, 0.317131, 0.310880,
2065 0.360456, 0.294052, 0.342694, 0.327166, 0.336797, 0.298040,
2066 0.295767, 0.260053, 0.325544, 0.335310, 0.320182, 0.301072,
2067 0.313117, 0.283407, 0.299206, 0.293525, 0.305067, 0.255978,
2068 0.327055, 0.316382, 0.317700, 0.278993, 0.283120, 0.314000,
2069 0.274396, 0.291208, 0.348813, 0.319603, 0.313076, 0.289155,
2070 0.343988, 0.311426, 0.322896, 0.328726, 0.337062, 0.389093,
2071 0.284122, 0.312184, 0.304008, 0.319170, 0.320778, 0.288390,
2072 0.337272, 0.356273, 0.343310, 0.312209, 0.330709, 0.297977,
2073 0.346146, 0.369162, 0.324385, 0.339831, 0.337037, 0.318739,
2074 0.343157, 0.277720, 0.368407, 0.321330, 0.338997, 0.314220,
2075 0.328861, 0.321824, 0.328013, 0.356925, 0.359144, 0.296314,
2076 0.345415, 0.396711, 0.347032, 0.294928, 0.343799, 0.322331,
2077 0.328656, 0.326098, 0.337338, 0.337038, 0.300179, 0.351391,
2078 0.324337, 0.330896, 0.302842, 0.310522, 0.337052, 0.359989,
2079 0.383250, 0.359355, 0.315382, 0.333113, 0.342598, 0.355348,
2080 0.320751, 0.320475, 0.351762, 0.351475, 0.338358, 0.326153,
2081 0.302507, 0.340048, 0.318685, 0.381646, 0.339320, 0.299453,
2082 0.426599, 0.393515, 0.353929, 0.328435, 0.413976, 0.292558,
2083 0.379340, 0.358344, 0.409259, 0.313821, 0.336675, 0.324521,
2084 0.408382, 0.346273, 0.312939, 0.362453, 0.343152, 0.330577,
2085 0.332831, 0.353299, 0.347745, 0.334818, 0.332234, 0.385585,
2086 0.395483, 0.395316, 0.326972, 0.349434, 0.317604, 0.328980,
2087 0.375056, 0.317290, 0.357083, 0.346165, 0.310444, 0.356873,
2088 0.359523, 0.365308, 0.365122, 0.383685, 0.370975, 0.396928,
2089 0.407654, 0.307755, 0.323033, 0.350580, 0.345231, 0.342462,
2090 0.400000, 0.318309, 0.403570, 0.322856, 0.383053, 0.422252,
2091 0.386112, 0.364314, 0.434375, 0.334629};
2092 wDD = 0.379611;
2093 wND = 0.496961;
2094 wSD = -1;
2095
2096     Mmin = bin[0];
2097     Mmax = bin[nbin];
2098     if(M<Mmin || M>Mmax) return kTRUE;
2099
2100     Int_t ibin=nbin-1;
2101     for(Int_t i=1; i<=nbin; i++)
2102       if(M<=bin[i]) {
2103         ibin=i-1;
2104         //      printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2105         break;
2106       }
2107     wSD=w[ibin];
2108     return kTRUE;
2109   }
2110  else if(TMath::Abs(fEnergyCMS-2760)<1 ){
2111
2112 const Int_t nbin=400;
2113 Double_t bin[]={
2114 1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
2115 4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
2116 7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
2117 10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
2118 13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
2119 15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
2120 18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
2121 21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
2122 24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
2123 27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
2124 30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
2125 33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
2126 36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
2127 39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
2128 42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
2129 45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
2130 48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
2131 51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
2132 54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
2133 57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
2134 60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
2135 63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
2136 66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
2137 69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
2138 72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
2139 75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
2140 78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
2141 81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
2142 84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
2143 87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
2144 90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
2145 93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
2146 96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
2147 99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2148 102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2149 105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2150 108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2151 111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2152 114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2153 117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2154 120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2155 123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2156 126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2157 129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2158 132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2159 135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2160 138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2161 141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2162 144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2163 147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2164 150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2165 153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2166 156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2167 159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2168 162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2169 165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2170 168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2171 171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2172 174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2173 177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2174 180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2175 183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2176 186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2177 189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2178 192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2179 195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2180 198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2181 Double_t w[]={
2182 1.000000, 0.692593, 0.673384, 0.666273, 0.657285, 0.637723,
2183 0.625881, 0.643590, 0.606100, 0.589007, 0.567824, 0.578705,
2184 0.538530, 0.517937, 0.528278, 0.515918, 0.539461, 0.466186,
2185 0.489869, 0.468402, 0.465017, 0.453336, 0.460769, 0.474638,
2186 0.456347, 0.434471, 0.427478, 0.435435, 0.410934, 0.366431,
2187 0.382240, 0.379513, 0.394249, 0.386837, 0.353103, 0.382138,
2188 0.377497, 0.389479, 0.378736, 0.347933, 0.354605, 0.352077,
2189 0.324443, 0.358792, 0.339968, 0.359052, 0.330734, 0.318291,
2190 0.333703, 0.358644, 0.335819, 0.332213, 0.309051, 0.309975,
2191 0.331626, 0.304407, 0.309819, 0.312097, 0.312462, 0.320411,
2192 0.280401, 0.302311, 0.315863, 0.281479, 0.310003, 0.296911,
2193 0.313676, 0.281071, 0.294163, 0.306500, 0.283462, 0.274867,
2194 0.307149, 0.270555, 0.282264, 0.287373, 0.307849, 0.278675,
2195 0.286990, 0.278269, 0.300105, 0.286799, 0.265674, 0.275140,
2196 0.285702, 0.257352, 0.267714, 0.248204, 0.252220, 0.255678,
2197 0.282946, 0.268257, 0.282375, 0.262675, 0.275564, 0.248345,
2198 0.236259, 0.291914, 0.259936, 0.241338, 0.267389, 0.285044,
2199 0.289419, 0.253594, 0.284568, 0.231840, 0.260008, 0.268527,
2200 0.275363, 0.224115, 0.281260, 0.257913, 0.295152, 0.264399,
2201 0.232287, 0.282533, 0.223431, 0.255756, 0.244471, 0.221695,
2202 0.272450, 0.284244, 0.253682, 0.270717, 0.275403, 0.240323,
2203 0.245081, 0.241859, 0.216340, 0.244789, 0.220291, 0.238478,
2204 0.224691, 0.244058, 0.266117, 0.271478, 0.242012, 0.267321,
2205 0.248494, 0.253343, 0.255606, 0.235458, 0.241079, 0.233223,
2206 0.226813, 0.259224, 0.234239, 0.258606, 0.210892, 0.238186,
2207 0.243271, 0.222678, 0.213437, 0.273939, 0.247966, 0.232548,
2208 0.263438, 0.222089, 0.272111, 0.248818, 0.244295, 0.238368,
2209 0.236908, 0.248776, 0.232604, 0.231194, 0.227117, 0.231152,
2210 0.282140, 0.229778, 0.232631, 0.261794, 0.216633, 0.253471,
2211 0.242157, 0.227406, 0.269335, 0.230547, 0.210618, 0.285872,
2212 0.248776, 0.229875, 0.242728, 0.227388, 0.220567, 0.222062,
2213 0.235950, 0.224087, 0.228895, 0.208287, 0.235999, 0.208696,
2214 0.230367, 0.267667, 0.220484, 0.233402, 0.233815, 0.250455,
2215 0.253120, 0.219556, 0.230980, 0.236661, 0.222395, 0.226111,
2216 0.198315, 0.210555, 0.202813, 0.208594, 0.235976, 0.221490,
2217 0.243059, 0.204901, 0.216987, 0.229039, 0.231466, 0.221975,
2218 0.231220, 0.253638, 0.250448, 0.260291, 0.328345, 0.205739,
2219 0.222014, 0.251513, 0.279427, 0.270506, 0.248409, 0.222472,
2220 0.291632, 0.227796, 0.248769, 0.276896, 0.214742, 0.200139,
2221 0.230693, 0.226031, 0.268900, 0.185160, 0.245353, 0.205843,
2222 0.231155, 0.219122, 0.214811, 0.199763, 0.274179, 0.217598,
2223 0.274988, 0.237244, 0.211820, 0.225459, 0.252799, 0.235948,
2224 0.224986, 0.245385, 0.237770, 0.213373, 0.229737, 0.215487,
2225 0.234453, 0.249684, 0.239435, 0.250422, 0.257194, 0.214762,
2226 0.212266, 0.228988, 0.253798, 0.201607, 0.239946, 0.205245,
2227 0.231670, 0.212774, 0.206768, 0.231563, 0.189388, 0.227926,
2228 0.227203, 0.237754, 0.221628, 0.211138, 0.203322, 0.200985,
2229 0.231780, 0.220294, 0.232686, 0.234243, 0.218264, 0.255870,
2230 0.223213, 0.238670, 0.213713, 0.213064, 0.246700, 0.233446,
2231 0.221503, 0.206767, 0.200722, 0.226179, 0.237425, 0.239229,
2232 0.238611, 0.240419, 0.247806, 0.215923, 0.205298, 0.232778,
2233 0.272312, 0.226773, 0.258103, 0.223287, 0.269404, 0.203398,
2234 0.223782, 0.204213, 0.229664, 0.234040, 0.228419, 0.203936,
2235 0.263686, 0.199141, 0.236127, 0.214058, 0.204611, 0.224324,
2236 0.292140, 0.190735, 0.235157, 0.213018, 0.257085, 0.190554,
2237 0.203197, 0.213044, 0.237023, 0.214243, 0.193562, 0.262403,
2238 0.206256, 0.221396, 0.233588, 0.256611, 0.249731, 0.226683,
2239 0.199330, 0.251026, 0.222596, 0.201941, 0.186374, 0.221038,
2240 0.196555, 0.222560, 0.299419, 0.231979, 0.242924, 0.198310,
2241 0.217628, 0.235458, 0.278595, 0.218624, 0.277305, 0.239109,
2242 0.205600, 0.253715, 0.221173, 0.218195, 0.277647, 0.241974,
2243 0.268748, 0.268128, 0.292258, 0.249420, 0.191034, 0.219506,
2244 0.216502, 0.250677, 0.193386, 0.201310, 0.259464, 0.255351,
2245 0.269628, 0.221063, 0.294079, 0.196726, 0.233634, 0.221870,
2246 0.216236, 0.197259, 0.247433, 0.272765, 0.294079, 0.236336,
2247 0.206396, 0.238524, 0.247846, 0.269519, 0.237141, 0.230611,
2248 0.201712, 0.242225, 0.255565, 0.258738};
2249 wDD = 0.512813;
2250 wND = 0.518820;
2251 wSD = -1;
2252
2253     Mmin = bin[0];
2254     Mmax = bin[nbin];
2255     if(M<Mmin || M>Mmax) return kTRUE;
2256
2257     Int_t ibin=nbin-1;
2258     for(Int_t i=1; i<=nbin; i++)
2259       if(M<=bin[i]) {
2260         ibin=i-1;
2261         //      printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2262         break;
2263       }
2264     wSD=w[ibin];
2265     return kTRUE;
2266   }
2267
2268
2269   else if(TMath::Abs(fEnergyCMS-7000)<1 ){
2270 const Int_t nbin=400;
2271 Double_t bin[]={
2272 1.080000, 1.577300, 2.074600, 2.571900, 3.069200, 3.566500,
2273 4.063800, 4.561100, 5.058400, 5.555700, 6.053000, 6.550300,
2274 7.047600, 7.544900, 8.042200, 8.539500, 9.036800, 9.534100,
2275 10.031400, 10.528700, 11.026000, 11.523300, 12.020600, 12.517900,
2276 13.015200, 13.512500, 14.009800, 14.507100, 15.004400, 15.501700,
2277 15.999000, 16.496300, 16.993600, 17.490900, 17.988200, 18.485500,
2278 18.982800, 19.480100, 19.977400, 20.474700, 20.972000, 21.469300,
2279 21.966600, 22.463900, 22.961200, 23.458500, 23.955800, 24.453100,
2280 24.950400, 25.447700, 25.945000, 26.442300, 26.939600, 27.436900,
2281 27.934200, 28.431500, 28.928800, 29.426100, 29.923400, 30.420700,
2282 30.918000, 31.415300, 31.912600, 32.409900, 32.907200, 33.404500,
2283 33.901800, 34.399100, 34.896400, 35.393700, 35.891000, 36.388300,
2284 36.885600, 37.382900, 37.880200, 38.377500, 38.874800, 39.372100,
2285 39.869400, 40.366700, 40.864000, 41.361300, 41.858600, 42.355900,
2286 42.853200, 43.350500, 43.847800, 44.345100, 44.842400, 45.339700,
2287 45.837000, 46.334300, 46.831600, 47.328900, 47.826200, 48.323500,
2288 48.820800, 49.318100, 49.815400, 50.312700, 50.810000, 51.307300,
2289 51.804600, 52.301900, 52.799200, 53.296500, 53.793800, 54.291100,
2290 54.788400, 55.285700, 55.783000, 56.280300, 56.777600, 57.274900,
2291 57.772200, 58.269500, 58.766800, 59.264100, 59.761400, 60.258700,
2292 60.756000, 61.253300, 61.750600, 62.247900, 62.745200, 63.242500,
2293 63.739800, 64.237100, 64.734400, 65.231700, 65.729000, 66.226300,
2294 66.723600, 67.220900, 67.718200, 68.215500, 68.712800, 69.210100,
2295 69.707400, 70.204700, 70.702000, 71.199300, 71.696600, 72.193900,
2296 72.691200, 73.188500, 73.685800, 74.183100, 74.680400, 75.177700,
2297 75.675000, 76.172300, 76.669600, 77.166900, 77.664200, 78.161500,
2298 78.658800, 79.156100, 79.653400, 80.150700, 80.648000, 81.145300,
2299 81.642600, 82.139900, 82.637200, 83.134500, 83.631800, 84.129100,
2300 84.626400, 85.123700, 85.621000, 86.118300, 86.615600, 87.112900,
2301 87.610200, 88.107500, 88.604800, 89.102100, 89.599400, 90.096700,
2302 90.594000, 91.091300, 91.588600, 92.085900, 92.583200, 93.080500,
2303 93.577800, 94.075100, 94.572400, 95.069700, 95.567000, 96.064300,
2304 96.561600, 97.058900, 97.556200, 98.053500, 98.550800, 99.048100,
2305 99.545400, 100.042700, 100.540000, 101.037300, 101.534600, 102.031900,
2306 102.529200, 103.026500, 103.523800, 104.021100, 104.518400, 105.015700,
2307 105.513000, 106.010300, 106.507600, 107.004900, 107.502200, 107.999500,
2308 108.496800, 108.994100, 109.491400, 109.988700, 110.486000, 110.983300,
2309 111.480600, 111.977900, 112.475200, 112.972500, 113.469800, 113.967100,
2310 114.464400, 114.961700, 115.459000, 115.956300, 116.453600, 116.950900,
2311 117.448200, 117.945500, 118.442800, 118.940100, 119.437400, 119.934700,
2312 120.432000, 120.929300, 121.426600, 121.923900, 122.421200, 122.918500,
2313 123.415800, 123.913100, 124.410400, 124.907700, 125.405000, 125.902300,
2314 126.399600, 126.896900, 127.394200, 127.891500, 128.388800, 128.886100,
2315 129.383400, 129.880700, 130.378000, 130.875300, 131.372600, 131.869900,
2316 132.367200, 132.864500, 133.361800, 133.859100, 134.356400, 134.853700,
2317 135.351000, 135.848300, 136.345600, 136.842900, 137.340200, 137.837500,
2318 138.334800, 138.832100, 139.329400, 139.826700, 140.324000, 140.821300,
2319 141.318600, 141.815900, 142.313200, 142.810500, 143.307800, 143.805100,
2320 144.302400, 144.799700, 145.297000, 145.794300, 146.291600, 146.788900,
2321 147.286200, 147.783500, 148.280800, 148.778100, 149.275400, 149.772700,
2322 150.270000, 150.767300, 151.264600, 151.761900, 152.259200, 152.756500,
2323 153.253800, 153.751100, 154.248400, 154.745700, 155.243000, 155.740300,
2324 156.237600, 156.734900, 157.232200, 157.729500, 158.226800, 158.724100,
2325 159.221400, 159.718700, 160.216000, 160.713300, 161.210600, 161.707900,
2326 162.205200, 162.702500, 163.199800, 163.697100, 164.194400, 164.691700,
2327 165.189000, 165.686300, 166.183600, 166.680900, 167.178200, 167.675500,
2328 168.172800, 168.670100, 169.167400, 169.664700, 170.162000, 170.659300,
2329 171.156600, 171.653900, 172.151200, 172.648500, 173.145800, 173.643100,
2330 174.140400, 174.637700, 175.135000, 175.632300, 176.129600, 176.626900,
2331 177.124200, 177.621500, 178.118800, 178.616100, 179.113400, 179.610700,
2332 180.108000, 180.605300, 181.102600, 181.599900, 182.097200, 182.594500,
2333 183.091800, 183.589100, 184.086400, 184.583700, 185.081000, 185.578300,
2334 186.075600, 186.572900, 187.070200, 187.567500, 188.064800, 188.562100,
2335 189.059400, 189.556700, 190.054000, 190.551300, 191.048600, 191.545900,
2336 192.043200, 192.540500, 193.037800, 193.535100, 194.032400, 194.529700,
2337 195.027000, 195.524300, 196.021600, 196.518900, 197.016200, 197.513500,
2338 198.010800, 198.508100, 199.005400, 199.502700, 200.000000};
2339 Double_t w[]={
2340 1.000000, 0.798640, 0.770197, 0.717990, 0.699434, 0.692093,
2341 0.620940, 0.655294, 0.636496, 0.633483, 0.600518, 0.594548,
2342 0.550759, 0.550413, 0.528926, 0.525877, 0.506701, 0.504117,
2343 0.486973, 0.492930, 0.461833, 0.488695, 0.438874, 0.460274,
2344 0.451551, 0.454424, 0.449270, 0.387571, 0.427556, 0.428740,
2345 0.401682, 0.402260, 0.408961, 0.395417, 0.387426, 0.378602,
2346 0.357826, 0.359125, 0.348494, 0.363710, 0.356117, 0.340363,
2347 0.337637, 0.396084, 0.341777, 0.340551, 0.348838, 0.344014,
2348 0.340468, 0.317654, 0.355584, 0.326023, 0.373416, 0.312298,
2349 0.326522, 0.290540, 0.335557, 0.318689, 0.327544, 0.319501,
2350 0.331754, 0.312728, 0.282263, 0.274937, 0.303867, 0.307820,
2351 0.289344, 0.268934, 0.288908, 0.290018, 0.291369, 0.295242,
2352 0.289067, 0.277685, 0.267957, 0.267559, 0.320229, 0.265060,
2353 0.305931, 0.305352, 0.262064, 0.281879, 0.287780, 0.270033,
2354 0.270814, 0.276667, 0.271531, 0.275283, 0.258189, 0.287969,
2355 0.251247, 0.301527, 0.267230, 0.245860, 0.293125, 0.253421,
2356 0.272396, 0.243637, 0.236206, 0.278452, 0.246544, 0.263165,
2357 0.230484, 0.231102, 0.258281, 0.244707, 0.270111, 0.248295,
2358 0.246942, 0.245592, 0.272766, 0.254243, 0.199647, 0.262590,
2359 0.226710, 0.243836, 0.214153, 0.233535, 0.235207, 0.234245,
2360 0.247698, 0.248379, 0.241463, 0.265230, 0.223242, 0.236191,
2361 0.252700, 0.231865, 0.228352, 0.229101, 0.237385, 0.246485,
2362 0.254706, 0.245565, 0.224841, 0.257301, 0.240968, 0.202379,
2363 0.236943, 0.241683, 0.220809, 0.219014, 0.213015, 0.244470,
2364 0.221042, 0.198996, 0.225295, 0.264466, 0.200600, 0.228143,
2365 0.250138, 0.216784, 0.268002, 0.275734, 0.218144, 0.229866,
2366 0.235443, 0.208909, 0.215067, 0.189743, 0.216741, 0.242686,
2367 0.200234, 0.218882, 0.245991, 0.222815, 0.219576, 0.209773,
2368 0.205247, 0.203855, 0.202534, 0.192536, 0.223447, 0.225810,
2369 0.220583, 0.248421, 0.223424, 0.206033, 0.203130, 0.221225,
2370 0.223763, 0.239216, 0.252311, 0.206893, 0.228461, 0.233591,
2371 0.201512, 0.242382, 0.215147, 0.232578, 0.207116, 0.239213,
2372 0.196215, 0.184367, 0.235135, 0.189768, 0.274084, 0.206267,
2373 0.203283, 0.248828, 0.250214, 0.217465, 0.232080, 0.215150,
2374 0.207936, 0.176789, 0.191338, 0.188655, 0.196181, 0.223691,
2375 0.254257, 0.216874, 0.213536, 0.221399, 0.192024, 0.212534,
2376 0.218169, 0.226635, 0.201191, 0.212700, 0.211634, 0.232353,
2377 0.223636, 0.208605, 0.249132, 0.183681, 0.221842, 0.187136,
2378 0.203772, 0.249575, 0.217713, 0.205193, 0.193941, 0.219707,
2379 0.264226, 0.182105, 0.207183, 0.220845, 0.236571, 0.182390,
2380 0.212914, 0.186266, 0.195361, 0.217665, 0.204527, 0.188804,
2381 0.222832, 0.191193, 0.201073, 0.185616, 0.195011, 0.200828,
2382 0.181434, 0.233707, 0.202925, 0.211992, 0.173100, 0.205258,
2383 0.182695, 0.216520, 0.218422, 0.209490, 0.211257, 0.215801,
2384 0.220652, 0.211409, 0.195731, 0.194957, 0.198888, 0.234926,
2385 0.221377, 0.229822, 0.176700, 0.172322, 0.212265, 0.206133,
2386 0.170355, 0.253305, 0.198688, 0.240043, 0.225384, 0.174729,
2387 0.195820, 0.200093, 0.196912, 0.212308, 0.200490, 0.240415,
2388 0.159744, 0.173686, 0.223853, 0.213604, 0.193779, 0.179328,
2389 0.180873, 0.237481, 0.179640, 0.235857, 0.202847, 0.167690,
2390 0.245093, 0.215504, 0.205848, 0.184408, 0.201626, 0.209651,
2391 0.236306, 0.217803, 0.188534, 0.187861, 0.161663, 0.221718,
2392 0.152044, 0.190412, 0.173505, 0.208995, 0.174575, 0.180734,
2393 0.247704, 0.203388, 0.194184, 0.169679, 0.182703, 0.213402,
2394 0.191808, 0.178604, 0.178116, 0.198452, 0.174687, 0.169809,
2395 0.222851, 0.156811, 0.229170, 0.210359, 0.178557, 0.248570,
2396 0.208536, 0.192571, 0.178912, 0.224505, 0.170822, 0.205492,
2397 0.330973, 0.160924, 0.203724, 0.213255, 0.205827, 0.187162,
2398 0.181252, 0.191723, 0.238106, 0.182398, 0.202358, 0.212066,
2399 0.201255, 0.168159, 0.185219, 0.229176, 0.158222, 0.164092,
2400 0.215405, 0.200724, 0.234811, 0.184488, 0.213112, 0.198577,
2401 0.219622, 0.160671, 0.179349, 0.206681, 0.206091, 0.165618,
2402 0.165203, 0.174442, 0.179287, 0.187318, 0.163481, 0.217440,
2403 0.160381, 0.177025, 0.204385, 0.163676, 0.210733, 0.186519,
2404 0.230701, 0.191764, 0.185119, 0.190770, 0.242987, 0.186812,
2405 0.202906, 0.161935, 0.182426, 0.197922, 0.181475, 0.155903,
2406 0.175006, 0.223482, 0.202706, 0.218108};
2407 wDD = 0.207705;
2408 wND = 0.289628;
2409 wSD = -1;
2410
2411     Mmin = bin[0];
2412     Mmax = bin[nbin];
2413
2414     if(M<Mmin || M>Mmax) return kTRUE;
2415
2416     Int_t ibin=nbin-1;
2417     for(Int_t i=1; i<=nbin; i++)
2418       if(M<=bin[i]) {
2419         ibin=i-1;
2420         //      printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
2421         break;
2422       }
2423     wSD=w[ibin];
2424     return kTRUE;
2425   }
2426
2427   return kFALSE;
2428 }
2429
2430 Bool_t AliGenPythia::IsFromHeavyFlavor(Int_t ipart)
2431 {
2432 // Check if this is a heavy flavor decay product
2433   TParticle *  part = (TParticle *) fParticles.At(ipart);
2434   Int_t mpdg = TMath::Abs(part->GetPdgCode());
2435   Int_t mfl  = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2436   //
2437   // Light hadron
2438   if (mfl >= 4 && mfl < 6) return kTRUE;
2439
2440   Int_t imo = part->GetFirstMother()-1;
2441   TParticle* pm = part;
2442   //
2443   // Heavy flavor hadron produced by generator
2444   while (imo >  -1) {
2445     pm  =  (TParticle*)fParticles.At(imo);
2446     mpdg = TMath::Abs(pm->GetPdgCode());
2447     mfl  = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
2448     if ((mfl > 3) && (mfl <6) && mpdg > 400) return kTRUE;
2449     imo = pm->GetFirstMother()-1;
2450   }
2451   return kFALSE;
2452 }
2453
2454 Bool_t AliGenPythia::CheckDetectorAcceptance(Float_t phi, Float_t eta, Int_t iparticle)
2455 {
2456   // check the eta/phi correspond to the detectors acceptance
2457   // iparticle is the index of the particle to be checked, for PHOS rotation case
2458   if     (fCheckPHOS   && IsInPHOS  (phi,eta,iparticle)) return kTRUE;
2459   else if(fCheckEMCAL  && IsInEMCAL (phi,eta)) return kTRUE;
2460   else if(fCheckBarrel && IsInBarrel(    eta)) return kTRUE;
2461   else                                         return kFALSE;
2462 }
2463
2464 Bool_t AliGenPythia::TriggerOnSelectedParticles(Int_t np)
2465 {
2466   // Select events with fragmentation photon, decay photon, pi0, eta or other hadrons going to PHOS or EMCAL or central barrel,
2467   // implemented primaryly for kPyJets, but extended to any kind of process.
2468
2469   //printf("Check: frag photon %d, pi0 %d, eta %d, electron %d, hadron %d, decay %d, PHOS %d, EMCAL %d, Barrel %d \n",
2470   //       fFragPhotonInCalo,fPi0InCalo, fEtaInCalo,fEleInCalo,fHadronInCalo,fDecayPhotonInCalo,fCheckPHOS,fCheckEMCAL, fCheckBarrel);
2471
2472   Bool_t ok = kFALSE;
2473   for (Int_t i=0; i< np; i++) {
2474
2475     TParticle* iparticle = (TParticle *) fParticles.At(i);
2476
2477     Int_t pdg          = iparticle->GetPdgCode();
2478     Int_t status       = iparticle->GetStatusCode();
2479     Int_t imother      = iparticle->GetFirstMother() - 1;
2480
2481     TParticle* pmother = 0x0;
2482     Int_t momStatus    = -1;
2483     Int_t momPdg       = -1;
2484     if(imother > 0 ){
2485       pmother = (TParticle *) fParticles.At(imother);
2486       momStatus    = pmother->GetStatusCode();
2487       momPdg       = pmother->GetPdgCode();
2488     }
2489
2490     ok = kFALSE;
2491
2492     //
2493     // Check the particle type: hadron (not pi0 or eta), electron, decay photon (from pi0 or eta or any), pi0 or eta
2494     //
2495     // Hadron
2496     if (fHadronInCalo && status == 1)
2497     {
2498       if(TMath::Abs(pdg) > 23 && pdg !=221 && pdg != 111) // avoid photons, electrons, muons, neutrinos and eta or pi0
2499         // (in case neutral mesons were declared stable)
2500         ok = kTRUE;
2501     }
2502     //Electron
2503     else if (fEleInCalo && status == 1 && TMath::Abs(pdg) == 11)
2504     {
2505         ok = kTRUE;
2506     }
2507     //Fragmentation photon
2508     else if (fFragPhotonInCalo && pdg == 22 && status == 1)
2509     {
2510       if(momStatus != 11) ok = kTRUE ;  // No photon from hadron decay
2511     }
2512     // Decay photon
2513     else if (fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay && pdg == 22) // pi0 status can be 1 or 11 depending on decay settings, work only for 11
2514     {
2515       if( momStatus == 11)
2516       {
2517         //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Decay photon! pdg %d, status %d, pt %2.2f, mom: pdg %d, pt %2.2f\n",
2518         //                                                   pdg,status,iparticle->Pt(),momPdg,pmother->Pt());
2519         ok = kTRUE ;  // photon from hadron decay
2520
2521         //In case only decays from pi0 or eta requested
2522         if(fPi0InCalo && momPdg!=111) ok = kFALSE;
2523         if(fEtaInCalo && momPdg!=221) ok = kFALSE;
2524       }
2525
2526     }
2527     // Pi0 or Eta particle
2528     else if ((fPi0InCalo || fEtaInCalo))
2529     {
2530       if(fDecayPhotonInCalo && !fForceNeutralMeson2PhotonDecay ) continue ;
2531
2532       if (fPi0InCalo && pdg == 111) // pi0 status can be 1 or 11 depending on decay settings
2533       {
2534         //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Pi0! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2535         ok = kTRUE;
2536       }
2537       else if (fEtaInCalo && pdg == 221)
2538       {
2539         //if(iparticle->Pt() > fTriggerParticleMinPt) printf("Eta! pdg %d, status %d, pt %2.2f\n",pdg,status,iparticle->Pt());
2540         ok = kTRUE;
2541       }
2542
2543     }// pi0 or eta
2544
2545     //
2546     // Check that the selected particle is in the calorimeter acceptance
2547     //
2548     if(ok && iparticle->Pt() > fTriggerParticleMinPt)
2549     {
2550       //Just check if the selected particle falls in the acceptance
2551       if(!fForceNeutralMeson2PhotonDecay )
2552       {
2553         //printf("\t Check acceptance! \n");
2554         Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2555         Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2556
2557         if(CheckDetectorAcceptance(phi,eta,i))
2558         {
2559           ok =kTRUE;
2560           AliDebug(1,Form("Selected trigger pdg %d, status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2561           //printf("\t Accept \n");
2562           break;
2563         }
2564         else ok = kFALSE;
2565       }
2566       //Mesons have several decay modes, select only those decaying into 2 photons
2567       else if(fForceNeutralMeson2PhotonDecay && (fPi0InCalo || fEtaInCalo))
2568       {
2569         // In case we want the pi0/eta trigger,
2570         // check the decay mode (2 photons)
2571
2572         //printf("\t Force decay 2 gamma\n");
2573
2574         Int_t ndaughters = iparticle->GetNDaughters();
2575         if(ndaughters != 2){
2576           ok=kFALSE;
2577           continue;
2578         }
2579
2580         TParticle*d1 = (TParticle *) fParticles.At(iparticle->GetDaughter(0)-1);
2581         TParticle*d2 = (TParticle *) fParticles.At(iparticle->GetDaughter(1)-1);
2582         if(!d1 || !d2) {
2583           ok=kFALSE;
2584           continue;
2585         }
2586
2587         //iparticle->Print();
2588         //d1->Print();
2589         //d2->Print();
2590
2591         Int_t pdgD1 = d1->GetPdgCode();
2592         Int_t pdgD2 = d2->GetPdgCode();
2593         //printf("\t \t 2 daughters pdg = %d - %d\n",pdgD1,pdgD2);
2594         //printf("mother %d - %d\n",d1->GetFirstMother(),d2->GetFirstMother());
2595
2596         if(pdgD1 != 22  || pdgD2 != 22){
2597           ok = kFALSE;
2598           continue;
2599         }
2600
2601         //printf("\t accept decay\n");
2602
2603         //Trigger on the meson, not on the daughter
2604         if(!fDecayPhotonInCalo){
2605
2606           Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
2607           Float_t eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax
2608
2609           if(CheckDetectorAcceptance(phi,eta,i))
2610           {
2611             //printf("\t Accept meson pdg %d\n",pdg);
2612             ok =kTRUE;
2613             AliDebug(1,Form("Selected trigger pdg %d (decay), status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2614             break;
2615           } else {
2616             ok=kFALSE;
2617             continue;
2618           }
2619         }
2620
2621         //printf("Check daughters acceptance\n");
2622
2623         //Trigger on the meson daughters
2624         //Photon 1
2625         Float_t phi = d1->Phi()*180./TMath::Pi(); //Convert to degrees
2626         Float_t eta =TMath::Abs(d1->Eta()); //in calos etamin=-etamax
2627         if(d1->Pt() > fTriggerParticleMinPt)
2628         {
2629           //printf("\t Check acceptance photon 1! \n");
2630           if(CheckDetectorAcceptance(phi,eta,i))
2631           {
2632             //printf("\t Accept Photon 1\n");
2633             ok =kTRUE;
2634             AliDebug(1,Form("Selected trigger pdg %d (decay), status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2635             break;
2636           }
2637           else ok = kFALSE;
2638         } // pt cut
2639         else  ok = kFALSE;
2640
2641         //Photon 2
2642         phi = d2->Phi()*180./TMath::Pi(); //Convert to degrees
2643         eta =TMath::Abs(d2->Eta()); //in calos etamin=-etamax
2644
2645         if(d2->Pt() > fTriggerParticleMinPt)
2646         {
2647           //printf("\t Check acceptance photon 2! \n");
2648           if(CheckDetectorAcceptance(phi,eta,i))
2649           {
2650             //printf("\t Accept Photon 2\n");
2651             ok =kTRUE;
2652             AliDebug(1,Form("Selected trigger pdg %d (decay), status %d, pt %2.2f, eta %2.2f, phi %2.2f\n",pdg,status,iparticle->Pt(), eta, phi));
2653             break;
2654           }
2655           else ok = kFALSE;
2656         } // pt cut
2657         else ok = kFALSE;
2658       } // force 2 photon daughters in pi0/eta decays
2659       else ok = kFALSE;
2660     } else ok = kFALSE; // check acceptance
2661   } // primary loop
2662
2663   //
2664   // If requested, rotate the particles event in phi to enhance/speed PHOS selection
2665   // A particle passing all trigger conditions except phi position in PHOS, is used as reference
2666   //
2667   if(fCheckPHOSeta)
2668   {
2669     RotatePhi(ok);
2670   }
2671
2672   return ok;
2673 }
2674
2675 void AliGenPythia::SetSeed(UInt_t seed)
2676 {
2677   GetRandom()->SetSeed(seed);
2678 }