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