PYTHIA6/AliGenPythia.cxx

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