PYTHIA6/AliGenPythia.cxx

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