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