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