Intial state radiation jets are also quenched.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 24 May 2004 07:05:57 +0000 (07:05 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 24 May 2004 07:05:57 +0000 (07:05 +0000)
PYTHIA6/AliPythia.cxx

index 03f885b137ce697931b6760f8bc3769887afe041..c0a9fa5a22dad964541dfba9a156091f4df20bea 100644 (file)
@@ -697,56 +697,81 @@ void  AliPythia::Quench()
     static Float_t eMean = 0.;
     static Int_t   icall = 0;
     
-    Double_t p0[2][5];
-    Double_t p1[2][5];
-    Double_t p2[2][5];
-    Int_t   klast[2] = {-1, -1};
+    Double_t p0[4][5];
+    Double_t p1[4][5];
+    Double_t p2[4][5];
+    Int_t   klast[4] = {-1, -1, -1, -1};
 
     Int_t numpart   = fPyjets->N;
     Double_t px = 0., py = 0., pz = 0., e = 0., m = 0., p = 0., pt = 0., theta = 0.;
-    Double_t pxq[2], pyq[2], pzq[2], eq[2], yq[2], mq[2], pq[2], phiq[2], thetaq[2], ptq[2];
-    Bool_t  quenched[2];
+    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 phi;
-    Double_t zInitial[2], wjtKick[2];
-    Int_t nGluon[2];
+    Double_t zInitial[4], wjtKick[4];
+    Int_t nGluon[4];
     
     Int_t   imo, kst, pdg;
 //
-//  Primary partons
+//  Sore information about Primary partons
+//
+//  j =
+//  0, 1 partons from hard scattering
+//  2, 3 partons from initial state radiation
+// 
+    for (Int_t i = 2; i <= 7; i++) {
+       Int_t j = 0;
+       // Skip gluons that participate in hard scattering
+       if (i == 4 || i == 5) continue;
+       // Gluons from hard Scattering
+       if (i == 6 || i == 7) {
+           j = i - 6;
+           pxq[j]    = fPyjets->P[0][i];
+           pyq[j]    = fPyjets->P[1][i];
+           pzq[j]    = fPyjets->P[2][i];
+           eq[j]     = fPyjets->P[3][i];
+           mq[j]     = fPyjets->P[4][i];
+       } else {
+           // Gluons from initial state radiation
+           //
+           // Obtain 4-momentum vector from difference between original parton and parton after gluon 
+           // radiation. Energy is calculated independently because initial state radition does not 
+           // conserve strictly momentum and energy for each partonic system independently.
+           //
+           // Not very clean. Should be improved !
+           //
+           //
+           j = i;
+           pxq[j]    = fPyjets->P[0][i] - fPyjets->P[0][i+2];
+           pyq[j]    = fPyjets->P[1][i] - fPyjets->P[1][i+2];
+           pzq[j]    = fPyjets->P[2][i] - fPyjets->P[2][i+2];
+           mq[j]     = fPyjets->P[4][i];
+           eq[j]     = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j] + mq[j] * mq[j]);
+       }
+//
+//  Calculate some kinematic variables
 //
-
-    
-    
-    for (Int_t i = 6; i <= 7; i++) {
-       Int_t j = i - 6;
-       
-       pxq[j]    = fPyjets->P[0][i];
-       pyq[j]    = fPyjets->P[1][i];
-       pzq[j]    = fPyjets->P[2][i];
-       eq[j]     = fPyjets->P[3][i];
-       mq[j]     = fPyjets->P[4][i];
        yq[j]     = 0.5 * TMath::Log((eq[j] + pzq[j] + 1.e-14) / (eq[j] - pzq[j] + 1.e-14));
        pq[j]     = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j] + pzq[j] * pzq[j]);
        phiq[j]   = TMath::Pi()+TMath::ATan2(-pyq[j], -pxq[j]);
        ptq[j]    = TMath::Sqrt(pxq[j] * pxq[j] + pyq[j] * pyq[j]);
        thetaq[j] = TMath::ATan2(ptq[j], pzq[j]);
        phi       = phiq[j];
-       
-       // Quench only central jets
-       if (TMath::Abs(yq[j]) > 2.5) {
+       //
+       // Quench only central jets and with E > 10.
+       //
+       if (TMath::Abs(yq[j]) > 2.5 || eq[j] < 10.) {
            zInitial[j] = 0.;
        } else {
            pdg =  fPyjets->K[1][i];
-           
            // Get length in nucleus
            Double_t l;
            fGlauber->GetLengthsForPythia(1, &phi, &l, -1.);
            //
-           // Energy loss for given length and parton typr 
+           // Energy loss for given length and parton type 
            Int_t itype = (pdg == 21) ? 2 : 1;
        
            Double_t eloss   = fQuenchingWeights->GetELossRandom(itype, l, eq[j]);
-           if (eq[j] > 80. && TMath::Abs(yq[j]) < 0.5) {
+           if (eq[j] > 40. && TMath::Abs(yq[j]) < 0.5) {
                icall ++;
                eMean += eloss;
            }
@@ -767,19 +792,19 @@ 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], phi, l, eloss, wjtKick[j], eMean / Float_t(icall+1), yq[j]);
            
-           zInitial[j] = 1.;
-           while (zInitial[j] >= 0.95)  zInitial[j] = gRandom->Exp(0.2);
+//         zInitial[j] = 0.8;
+//         while (zInitial[j] >= 0.95)  zInitial[j] = gRandom->Exp(0.2);
        }
        
        quenched[j] = (zInitial[j] > 0.01);
     } // primary partons
-  
+    
     Double_t pNew[1000][4];
     Int_t    kNew[1000];
     Int_t icount = 0;
 //
 //  System Loop    
-    for (Int_t isys = 0; isys < 2; isys++) {
+    for (Int_t isys = 0; isys < 4; isys++) {
 //      Skip to next system if not quenched.
        if (!quenched[isys]) continue;
        
@@ -798,9 +823,6 @@ void  AliPythia::Quench()
 // Loop on radiation events
 
        for (Int_t iglu = 0; iglu < nGluon[isys]; iglu++) {
-           Double_t zHeavy = zInitial[isys];
-//
-
            while (1) {
                icount = 0;
                for (Int_t k = 0; k < 4; k++)
@@ -821,8 +843,12 @@ void  AliPythia::Quench()
 //      Quarks and gluons only
                    if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
 //      Particles from hard scattering only
+                   
                    if (imo > 8 && imo < 1000) imo = fPyjets->K[2][imo - 1];
-                   if (imo != (isys + 7) && (imo % 1000)  != (isys + 7)) continue;
+                   Int_t imom = imo % 1000;
+                   if ((isys == 0 || isys == 1) && ((imom != (isys + 7)))) continue;
+                   if ((isys == 2 || isys == 3) && ((imom != (isys + 1)))) continue;               
+                   
                    
 //      Skip comment lines
                    if (kst != 1 && kst != 2) continue;
@@ -839,9 +865,8 @@ void  AliPythia::Quench()
                    theta = TMath::ATan2(pt, pz);
                
 //
-//      Save 4-momentum sum for balancing      
-                   Int_t index = imo - 7;
-                   if (index >=  1000) index = imo % 1000 - 7;
+//      Save 4-momentum sum for balancing
+                   Int_t index = isys;
                    
                    p0[index][0] += px;
                    p0[index][1] += py;
@@ -854,18 +879,19 @@ void  AliPythia::Quench()
 //      Fractional energy loss
                    Double_t z = zInitial[index];
                    
+                   
 //      Don't fully quench radiated gluons
 //
                    if (imo > 1000) {
 //      This small factor makes sure that the gluons are not too close in phase space to avoid recombination
 //
 
-                       z = 0.05;
+                       z = 0.02;
                    }
+//                 printf("z: %d %f\n", imo, z);
+                   
 
 //
-
-                   if (m > 0.) z = zHeavy;
                    
                    //
                    //