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