Store and communicate event specific quenching parameters.

diff --git a/PYTHIA6/AliPythia.cxx b/PYTHIA6/AliPythia.cxx
index 75e8d30105d08c4b5edb3454051b3d7bd6ed985e..43ab97e06d8af9e7a3ec3a34751cfb2baf379881 100644 (file)
--- a/PYTHIA6/AliPythia.cxx
+++ b/PYTHIA6/AliPythia.cxx
@@ -714,10 +714,11 @@ void  AliPythia::Quench()
     Double_t px = 0., py = 0., pz = 0., e = 0., m = 0., p = 0., pt = 0., theta = 0., phi = 0.;
     Double_t pxq[4], pyq[4], pzq[4], eq[4], yq[4], mq[4], pq[4], phiq[4], thetaq[4], ptq[4];
     Bool_t  quenched[4];
-    Double_t zInitial[4], wjtKick[4];
+    Double_t wjtKick[4];
     Int_t nGluon[4];
     Int_t qPdg[4];
     Int_t   imo, kst, pdg;
+    
 //
 //  Sore information about Primary partons
 //
@@ -767,9 +768,9 @@ void  AliPythia::Quench()
   
     Double_t int0[4];
     Double_t int1[4];
-
-    fGlauber->GetI0I1ForPythia(4, phiq, int0, int1, 15.);
     
+    fGlauber->GetI0I1ForPythiaAndXY(4, phiq, int0, int1, fXJet, fYJet, 15.);
+
     for (Int_t j = 0; j < 4; j++) {
        //
        // Quench only central jets and with E > 10.
@@ -780,7 +781,7 @@ void  AliPythia::Quench()
        Double_t eloss = fQuenchingWeights->GetELossRandomKFast(itype, int0[j], int1[j], eq[j]);
 
        if (TMath::Abs(yq[j]) > 2.5 || eq[j] < 10.) {
-           zInitial[j] = 0.;
+           fZQuench[j] = 0.;
        } else {
            if (eq[j] > 40. && TMath::Abs(yq[j]) < 0.5) {
                icall ++;
@@ -792,11 +793,11 @@ void  AliPythia::Quench()
            wjtKick[j] = TMath::Sqrt(l *  fQuenchingWeights->CalcQk(int0[j], int1[j]));
            //
            // Fractional energy loss
-           zInitial[j] = eloss / eq[j];
+           fZQuench[j] = eloss / eq[j];
            //
            // Avoid complete loss
            //
-           if (zInitial[j] == 1.) zInitial[j] = 0.95;
+           if (fZQuench[j] == 1.) fZQuench[j] = 0.95;
            //
            // Some debug printing
 
@@ -804,25 +805,29 @@ void  AliPythia::Quench()
 //         printf("Initial parton # %3d, Type %3d Energy %10.3f Phi %10.3f Length %10.3f Loss %10.3f Kick %10.3f Mean: %10.3f %10.3f\n", 
 //                j, itype, eq[j], phiq[j], l, eloss, wjtKick[j], eMean / Float_t(icall+1), yq[j]);
            
-//         zInitial[j] = 0.8;
-//         while (zInitial[j] >= 0.95)  zInitial[j] = gRandom->Exp(0.2);
+//         fZQuench[j] = 0.8;
+//         while (fZQuench[j] >= 0.95)  fZQuench[j] = gRandom->Exp(0.2);
        }
        
-       quenched[j] = (zInitial[j] > 0.01);
+       quenched[j] = (fZQuench[j] > 0.01);
     } // primary partons
     
+    
+
     Double_t pNew[1000][4];
     Int_t    kNew[1000];
     Int_t icount = 0;
+    Double_t zquench[4];
+    
 //
 //  System Loop    
     for (Int_t isys = 0; isys < 4; isys++) {
 //      Skip to next system if not quenched.
        if (!quenched[isys]) continue;
        
-       nGluon[isys]   = 1 + Int_t(zInitial[isys] / (1. - zInitial[isys]));
+       nGluon[isys]   = 1 + Int_t(fZQuench[isys] / (1. - fZQuench[isys]));
        if (nGluon[isys] > 6) nGluon[isys] = 6;
-       zInitial[isys] = 1. - TMath::Power(1. - zInitial[isys], 1./Double_t(nGluon[isys]));
+       zquench[isys] = 1. - TMath::Power(1. - fZQuench[isys], 1./Double_t(nGluon[isys]));
        wjtKick[isys]  = wjtKick[isys] / TMath::Sqrt(Double_t(nGluon[isys]));
 
 
@@ -889,7 +894,7 @@ void  AliPythia::Quench()
                    
 //
 //      Fractional energy loss
-                   Double_t z = zInitial[index];
+                   Double_t z = zquench[index];
                    
                    
 //      Don't fully quench radiated gluons
@@ -1004,7 +1009,7 @@ void  AliPythia::Quench()
                    p2[isys][4] = TMath::Sqrt(p2[isys][4]);
                    break;
                } else {
-                   printf("Warning negative mass squared in system %d %f ! \n", isys, zInitial[isys]);
+                   printf("Warning negative mass squared in system %d %f ! \n", isys, zquench[isys]);
                    printf("4-Momentum: %10.3e %10.3e %10.3e %10.3e %10.3e \n", p2[isys][0], p2[isys][1], p2[isys][2], p2[isys][3], p2[isys][4]);
                    if (p2[isys][4] < -0.01) {
                        printf("Negative mass squared !\n");
@@ -1208,3 +1213,13 @@ void AliPythia::Pyquen(Double_t a, Int_t ibf, Double_t b)
     // Igor Lokthine's quenching routine
     pyquen(a, ibf, b);
 }
+
+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];
+
+}
+

diff --git a/PYTHIA6/AliPythia.h b/PYTHIA6/AliPythia.h
index 8ad7ae5235e66a8e1a5d96c7c0c497cb0bf45267..03aae7f95ddf04f84a99578f9b1f89b3539bc04c 100644 (file)
--- a/PYTHIA6/AliPythia.h
+++ b/PYTHIA6/AliPythia.h
@@ -58,6 +58,8 @@ class AliPythia : public TPythia6, public AliRndm
     virtual void Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez);
     virtual void InitQuenching(Float_t bmin, Float_t bmax, Float_t k, Int_t iECMethod);
     virtual void Pyquen(Double_t a, Int_t ibf, Double_t b);
+    virtual void GetQuenchingParameters(Double_t& xp, Double_t& yp, Double_t z[4]);
+    
     
     // return instance of the singleton
     static  AliPythia* Instance();
@@ -68,6 +70,9 @@ class AliPythia : public TPythia6, public AliRndm
     StrucFunc_t           fStrucFunc;         // Structure function
     Int_t                 fDefMDCY[501];      //  ! Default decay switches per particle
     Int_t                 fDefMDME[2001];     //  ! Default decay switches per mode
+    Double_t              fZQuench[4];        //  ! Quenching fractions for this even
+    Double_t              fXJet;              //  ! Jet production point X
+    Double_t              fYJet;              //  ! Jet production point Y
     AliFastGlauber*       fGlauber;           //  ! The Glauber model
     AliQuenchingWeights*  fQuenchingWeights;  //  ! The Quenching Weights model
     static AliPythia*     fgAliPythia;        // Pointer to single instance