PYTHIA6/AliGenPythia.cxx

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