Store and communicate event specific quenching parameters.
[u/mrichter/AliRoot.git] / PYTHIA6 / AliPythia.cxx
index 75e8d30105d08c4b5edb3454051b3d7bd6ed985e..43ab97e06d8af9e7a3ec3a34751cfb2baf379881 100644 (file)
@@ -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];
+
+}
+