PYTHIA6/AliGenPythia.cxx

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