PYTHIA6/AliGenPythia.cxx

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