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