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