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