PYTHIA6/AliGenPythia.cxx

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