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