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