PYTHIA6/AliGenPythia.cxx

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