+
+void AliPythia::Pyquen(Double_t a, Int_t ibf, Double_t b)
+{
+ // Igor Lokthine's quenching routine
+ // http://lokhtin.web.cern.ch/lokhtin/pyquen/pyquen.txt
+
+ pyquen(a, ibf, b);
+}
+
+void AliPythia::SetPyquenParameters(Double_t t0, Double_t tau0, Int_t nf, Int_t iengl, Int_t iangl)
+{
+ // Set the parameters for the PYQUEN package.
+ // See comments in PyquenCommon.h
+
+
+ PYQPAR.t0 = t0;
+ PYQPAR.tau0 = tau0;
+ PYQPAR.nf = nf;
+ PYQPAR.iengl = iengl;
+ PYQPAR.iangl = iangl;
+}
+
+
+void AliPythia::Pyevnw()
+{
+ // New multiple interaction scenario
+ pyevnw();
+}
+
+void AliPythia::Pyshowq(Int_t ip1, Int_t ip2, Double_t qmax)
+{
+ // Call medium-modified Pythia jet reconstruction algorithm
+ //
+ pyshowq(ip1, ip2, qmax);
+}
+ void AliPythia::Qpygin0()
+ {
+ // New multiple interaction scenario
+ qpygin0();
+ }
+
+void AliPythia::GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4])
+{
+ // Return event specific quenching parameters
+ xp = fXJet;
+ yp = fYJet;
+ for (Int_t i = 0; i < 4; i++) z[i] = fZQuench[i];
+
+}
+
+void AliPythia::ConfigHeavyFlavor()
+{
+ //
+ // Default configuration for Heavy Flavor production
+ //
+ // All QCD processes
+ //
+ SetMSEL(1);
+
+
+ if (fItune < 0) {
+ // No multiple interactions
+ SetMSTP(81,0);
+ SetPARP(81, 0.);
+ SetPARP(82, 0.);
+ }
+ // Initial/final parton shower on (Pythia default)
+ SetMSTP(61,1);
+ SetMSTP(71,1);
+
+ // 2nd order alpha_s
+ SetMSTP(2,2);
+
+ // QCD scales
+ SetMSTP(32,2);
+ SetPARP(34,1.0);
+}
+
+void AliPythia::AtlasTuning()
+{
+ //
+ // Configuration for the ATLAS tuning
+ if (fItune > -1) return;
+ printf("ATLAS TUNE \n");
+
+ SetMSTP(51, AliStructFuncType::PDFsetIndex(kCTEQ5L)); // CTEQ5L pdf
+ SetMSTP(81,1); // Multiple Interactions ON
+ SetMSTP(82,4); // Double Gaussian Model
+ SetPARP(81,1.9); // Min. pt for multiple interactions (default in 6.2-14)
+ SetPARP(82,1.8); // [GeV] PT_min at Ref. energy
+ SetPARP(89,1000.); // [GeV] Ref. energy
+ SetPARP(90,0.16); // 2*epsilon (exponent in power law)
+ SetPARP(83,0.5); // Core density in proton matter distribution (def.value)
+ SetPARP(84,0.5); // Core radius
+ SetPARP(85,0.33); // Regulates gluon prod. mechanism
+ SetPARP(86,0.66); // Regulates gluon prod. mechanism
+ SetPARP(67,1); // Regulates Initial State Radiation
+}
+
+void AliPythia::AtlasTuningMC09()
+{
+ //
+ // Configuration for the ATLAS tuning
+ if (fItune > -1) return;
+ printf("ATLAS New TUNE MC09\n");
+ SetMSTP(81,21); // treatment for MI, ISR, FSR and beam remnants: MI on, new model
+ SetMSTP(82, 4); // Double Gaussian Model
+ SetMSTP(52, 2); // External PDF
+ SetMSTP(51, 20650); // MRST LO*
+
+
+ SetMSTP(70, 0); // (was 2: def manual 1, def code 0) virtuality scale for ISR
+ SetMSTP(72, 1); // (was 0: def 1) maximum scale for FSR
+ SetMSTP(88, 1); // (was 0: def 1) strategy for qq junction to di-quark or baryon in beam remnant
+ SetMSTP(90, 0); // (was 1: def 0) strategy of compensate the primordial kT
+
+ SetPARP(78, 0.3); // the amount of color reconnection in the final state
+ SetPARP(80, 0.1); // probability of color partons kicked out from beam remnant
+ SetPARP(82, 2.3); // [GeV] PT_min at Ref. energy
+ SetPARP(83, 0.8); // Core density in proton matter distribution (def.value)
+ SetPARP(84, 0.7); // Core radius
+ SetPARP(90, 0.25); // 2*epsilon (exponent in power law)
+ SetPARJ(81, 0.29); // (was 0.14: def 0.29) Labmda value in running alpha_s for parton showers
+
+ SetMSTP(95, 6);
+ SetPARJ(41, 0.3); // a and b parameters of the symmm. Lund FF
+ SetPARJ(42, 0.58);
+ SetPARJ(46, 0.75); // mod. of the Lund FF for heavy end-point quarks
+ SetPARP(89,1800.); // [GeV] Ref. energy
+}
+
+AliPythia& AliPythia::operator=(const AliPythia& rhs)
+{
+// Assignment operator
+ rhs.Copy(*this);
+ return *this;
+}
+
+ void AliPythia::Copy(TObject&) const
+{
+ //
+ // Copy
+ //
+ Fatal("Copy","Not implemented!\n");
+}
+
+void AliPythia::DalitzDecays()
+{
+
+ //
+ // Replace all omega dalitz decays with the correct matrix element decays
+ //
+ Int_t nt = fPyjets->N;
+ for (Int_t i = 0; i < nt; i++) {
+ if (fPyjets->K[1][i] != 223) continue;
+ Int_t fd = fPyjets->K[3][i] - 1;
+ Int_t ld = fPyjets->K[4][i] - 1;
+ if (fd < 0) continue;
+ if ((ld - fd) != 2) continue;
+ if ((fPyjets->K[1][fd] != 111) ||
+ ((TMath::Abs(fPyjets->K[1][fd+1]) != 11) && (TMath::Abs(fPyjets->K[1][fd+1]) != 13)))
+ continue;
+ TLorentzVector omega(fPyjets->P[0][i], fPyjets->P[1][i], fPyjets->P[2][i], fPyjets->P[3][i]);
+ Int_t pdg = TMath::Abs(fPyjets->K[1][fd+1]);
+ fOmegaDalitz.Decay(pdg, &omega);
+ for (Int_t j = 0; j < 3; j++) {
+ for (Int_t k = 0; k < 4; k++) {
+ TLorentzVector vec = (fOmegaDalitz.Products())[2-j];
+ fPyjets->P[k][fd+j] = vec[k];
+ }
+ }
+ }
+}