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