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