PYTHIA6/AliGenPythia.cxx

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