Independent fragmentation of radiated gluons.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 22 Apr 2004 09:22:43 +0000 (09:22 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 22 Apr 2004 09:22:43 +0000 (09:22 +0000)
PYTHIA6/AliPythia.cxx

index 85fd6d166ca84df47da38eee4665a75a750d509f..a39e3ed85964e8cd42df9aaf1ec56b9183812f73 100644 (file)
@@ -531,17 +531,6 @@ void  AliPythia::Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)
 {
 //  Call Pythia jet reconstruction algorithm
 //
-    Int_t numpart   = fPyjets->N;
-    for (Int_t i = 0; i < numpart; i++)
-    {
-       if (fPyjets->K[2][i] == 7) ip1 = i+1;
-       if (fPyjets->K[2][i] == 8) ip2 = i+1;
-    }
-    
-
-    qmax = 2. * GetVINT(51);
-    printf("Pyshow %d %d %f", ip1, ip2, qmax);
-    
     pyshow(ip1, ip2, qmax);
 }
 
@@ -593,7 +582,8 @@ void  AliPythia::Quench()
 //
 //
 // 
-    const Int_t kGluons = 1;
+    static Float_t eMean = 0.;
+    static Int_t   icall = 0;
     
     Double_t p0[2][5];
     Double_t p1[2][5];
@@ -607,10 +597,14 @@ void  AliPythia::Quench()
     Bool_t  quenched[2];
     Double_t phi;
     Double_t zInitial[2], wjtKick[2];
+    Int_t nGluon[2];
+    
     Int_t   imo, kst, pdg;
 //
 //  Primary partons
 //
+
+    
     
     for (Int_t i = 6; i <= 7; i++) {
        Int_t j = i - 6;
@@ -620,7 +614,7 @@ void  AliPythia::Quench()
        pzq[j]    = fPyjets->P[2][i];
        eq[j]     = fPyjets->P[3][i];
        mq[j]     = fPyjets->P[4][i];
-       yq[j]     = 0.5 * TMath::Log((e + pz + 1.e-14) / (e - pz + 1.e-14));
+       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]);
@@ -639,9 +633,16 @@ void  AliPythia::Quench()
            //
            // Energy loss for given length and parton typr 
            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) {
+               icall ++;
+               eMean += eloss;
+           }
+           
            //
            // Extra pt
+           
            wjtKick[j] = TMath::Sqrt(l *  fQuenchingWeights->GetQTransport());
            //
            // Fractional energy loss
@@ -652,43 +653,43 @@ void  AliPythia::Quench()
            if (zInitial[j] == 1.) zInitial[j] = 0.95;
            //
            // Some debug printing
-           printf("Initial parton # %3d, Type %3d Energy %10.3f Phi %10.3f Length %10.3f Loss %10.3f Kick %10.3f\n", 
-                  j, itype, eq[j], phi, l, eloss, wjtKick[j]);
+           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);
        }
+       
        quenched[j] = (zInitial[j] > 0.01);
-    }
+    } // primary partons
   
-//
-// Radiated partons
-//
-    zInitial[0] = 1. - TMath::Power(1. - zInitial[0], 1./Double_t(kGluons));
-    zInitial[1] = 1. - TMath::Power(1. - zInitial[1], 1./Double_t(kGluons));
-    wjtKick[0]  = wjtKick[0] / TMath::Sqrt(Double_t(kGluons));
-    wjtKick[1]  = wjtKick[1] / TMath::Sqrt(Double_t(kGluons));
-//    this->Pylist(1);
-    
-//  Arrays to store particle 4-momenta to be changed
-//
-/*    
-    Double_t** pNew =  new Double_t* [numpart];
-    for (Int_t i = 0; i < numpart; i++) pNew[i] = new Double_t [4];
-    Int_t* kNew    = new Int_t [numpart];
-*/
-
     Double_t pNew[1000][4];
     Int_t    kNew[1000];
     Int_t icount = 0;
 //
-//  Radiation Loop    
-    for (Int_t iglu = 0; iglu < kGluons; iglu++) {
-       for (Int_t isys = 0; isys < 2; isys++) {
+//  System Loop    
+    for (Int_t isys = 0; isys < 2; isys++) {
 //      Skip to next system if not quenched.
+       if (!quenched[isys]) continue;
+       
+       nGluon[isys]   = 1 + Int_t(zInitial[isys] / (1. - zInitial[isys]));
+       if (nGluon[isys] > 6) nGluon[isys] = 6;
+       zInitial[isys] = 1. - TMath::Power(1. - zInitial[isys], 1./Double_t(nGluon[isys]));
+       wjtKick[isys]  = wjtKick[isys] / TMath::Sqrt(Double_t(nGluon[isys]));
 
-           
+
+       
+       Int_t igMin = -1;
+       Int_t igMax = -1;
+       Double_t pg[4] = {0., 0., 0., 0.};
+       
+//
+// Loop on radiation events
+
+       for (Int_t iglu = 0; iglu < nGluon[isys]; iglu++) {
            Double_t zHeavy = zInitial[isys];
-           
-           if (!quenched[isys]) continue;
 //
+
            while (1) {
                icount = 0;
                for (Int_t k = 0; k < 4; k++)
@@ -710,7 +711,7 @@ void  AliPythia::Quench()
                    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;
+                   if (imo != (isys + 7) && (imo % 1000)  != (isys + 7)) continue;
                    
 //      Skip comment lines
                    if (kst != 1 && kst != 2) continue;
@@ -729,30 +730,30 @@ void  AliPythia::Quench()
 //
 //      Save 4-momentum sum for balancing      
                    Int_t index = imo - 7;
-                   if (index >=  1000) index -= 1000;
+                   if (index >=  1000) index = imo % 1000 - 7;
                    
                    p0[index][0] += px;
                    p0[index][1] += py;
                    p0[index][2] += pz;
                    p0[index][3] += e;
-                   
-//      Don't quench radiated gluons
-//
-                   if (imo == 1000 + isys + 7) {
-                       p1[index][0] += px;
-                       p1[index][1] += py;
-                       p1[index][2] += pz;
-                       p1[index][3] += e;      
-                       continue;
-                   }
-                   
-//
                
                    klast[index] = i;
                    
 //
 //      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;
+                   }
+
+//
+
                    if (m > 0.) z = zHeavy;
                    
                    //
@@ -769,27 +770,7 @@ void  AliPythia::Quench()
                    //
                    // Kinematic limit on z
                    //
-                   
-                   if (m > 0.) {
-                       zmax = 1. - m / TMath::Sqrt(m * m + jt * jt);
-                       if (z > zmax) {
-                           printf("We have to put z to the kinematic limit %e %e \n", z, zmax);
-                           z = 0.9999 * zmax;
-                       } // z > zmax
-                       
-                       if (z < 0.01) {
-//
-//           If z is too small, there is no phase space for quenching
-//
-                           printf("No phase space for quenching ! %e  \n", z);
-                           p1[index][0] += px;
-                           p1[index][1] += py;
-                           p1[index][2] += pz;
-                           p1[index][3] += e;
-                           continue;
-                       }
-                   } // massive particles
-                   
+                   if (m > 0.) zmax = 1. - m / TMath::Sqrt(m * m + jt * jt);
                    //
                    // Change light-cone kinematics rel. to initial parton
                    //  
@@ -806,18 +787,24 @@ void  AliPythia::Quench()
                    // if mt very small (or sometimes even < 0 for numerical reasons) set it to 0
                    Double_t mt2New = eppzNew * empzNew;
                    if (mt2New < 1.e-8) mt2New = 0.;
-                   
-                   if (m * m > mt2New) {
-                       //
-                       // This should not happen 
-                       //
-                       Fatal("Quench()", "This should never happen %e %e %e!", m, eppzNew, empzNew);
-                       jtNew = 0;
+                   if (z < zmax) {
+                       if (m * m > mt2New) {
+                           //
+                           // This should not happen 
+                           //
+                           Fatal("Quench()", "This should never happen %e %e %e!", m, eppzNew, empzNew);
+                           jtNew = 0;
+                       } else {
+                           jtNew    = TMath::Sqrt(mt2New - m * m);
+                       }
                    } else {
-                       jtNew    = TMath::Sqrt(mt2New - m * m);
+                       // If pT is to small (probably a leading massive particle) we scale only the energy
+                       // This can cause negative masses of the radiated gluon
+                       // Let's hope for the best ...
+                       jtNew = jt;
+                       eNew  = TMath::Sqrt(plNew * plNew + mt2);
+                       
                    }
-               
-           
                    //
                    //     Calculate new px, py
                    //
@@ -864,42 +851,25 @@ void  AliPythia::Quench()
                    p2[isys][j] = p0[isys][j] - p1[isys][j];
                }
                p2[isys][4] = p2[isys][3] * p2[isys][3] - p2[isys][0] * p2[isys][0] - p2[isys][1] * p2[isys][1] - p2[isys][2] * p2[isys][2];
-               
                if (p2[isys][4] > 0.) {
                    p2[isys][4] = TMath::Sqrt(p2[isys][4]);
                    break;
                } else {
                    printf("Warning negative mass squared in system %d %f ! \n", isys, zInitial[isys]);
-                   
-                   printf("Kinematics %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]);
+                   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 ! Let's try to fix this by decreasing z\n");
-//                     this->Pylist(1);
-                       
+                       printf("Negative mass squared !\n");
+                       // Here we have to put the gluon back to mass shell
+                       // This will lead to a small energy imbalance
+                       p2[isys][4]  = 0.;
+                       p2[isys][3]  = TMath::Sqrt(p2[isys][0] * p2[isys][0] + p2[isys][1] * p2[isys][1] + p2[isys][2] * p2[isys][2]);
+                       break;
                    } else {
                        p2[isys][4] = 0.;
                        break;
                    }
                }
-               //
-               // jt-kick
-               //
                /*
-                 TVector3 v(p2[k][0], p2[k][1], p2[k][2]);
-                 v.RotateZ(-phiq[k]);
-                 v.RotateY(-thetaq[k]);
-                 Double_t px = v.X(); Double_t py = v.Y(); Double_t pz  = v.Z();          
-                 Double_t r       = AliPythiaRndm::GetPythiaRandom()->Rndm();
-                 Double_t  jtKick  = wjtKick[k] * TMath::Sqrt(-TMath::Log(r));
-                 Double_t phiKick = 2. * TMath::Pi() * AliPythiaRndm::GetPythiaRandom()->Rndm();
-                 px += jtKick * TMath::Cos(phiKick);
-                 py += jtKick * TMath::Sin(phiKick);
-                 TVector3 w(px, py, pz);
-                 w.RotateY(thetaq[k]);
-                 w.RotateZ(phiq[k]);
-                 p2[k][0] = w.X(); p2[k][1] = w.Y(); p2[k][2] = w.Z();
-                 p2[k][3] = TMath::Sqrt(p2[k][0] * p2[k][0] + p2[k][1] * p2[k][1] + p2[k][2] * p2[k][2] + p2[k][4] * p2[k][4]);
-               */
                zHeavy *= 0.98;
                printf("zHeavy lowered to %f\n", zHeavy);
                if (zHeavy < 0.01) {
@@ -908,8 +878,9 @@ void  AliPythia::Quench()
                    quenched[isys] = kFALSE;
                    break;
                }
-           } // iteration on z
-           
+               */
+           } // iteration on z (while)
+           
 //         Update  event record
            for (Int_t k = 0; k < icount; k++) {
 //             printf("%6d %6d %10.3e %10.3e %10.3e %10.3e\n", k, kNew[k], pNew[k][0],pNew[k][1], pNew[k][2], pNew[k][3] );
@@ -918,78 +889,62 @@ void  AliPythia::Quench()
                fPyjets->P[2][kNew[k]] = pNew[k][2];
                fPyjets->P[3][kNew[k]] = pNew[k][3];
            }
-       } // System         
-       //
-       // Add the gluons
-       //
-       Int_t ish = 0;    
-       for (Int_t i = 0; i < 2; i++) {
-           Int_t jmin, jmax, iGlu, iNew;
-           if (!quenched[i]) continue;
+           //
+           // Add the gluons
+           //
+           Int_t ish = 0;    
+           Int_t iGlu, jmin, jmax, iNew;
+           if (!quenched[isys]) continue;
 //
 //      Last parton from shower i
-           Int_t in = klast[i];
+           Int_t in = klast[isys];
 //
 //      Continue if no parton in shower i selected
            if (in == -1) continue;
 //  
 //      If this is the second initial parton and it is behind the first move pointer by previous ish
-           if (i == 1 && klast[1] > klast[0]) in += ish;
+           if (isys == 1 && klast[1] > klast[0]) in += ish;
 //
 //      Starting index
            
-           jmin = in - 1;
+//         jmin = in - 1;
 // How many additional gluons will be generated
            ish  = 1;
-           if (p2[i][4] > 0.05) ish = 2;
+           if (p2[isys][4] > 0.05) ish = 2;
 //
 //      Position of gluons
-           iGlu = in;
-           iNew = in + ish;
-           jmax = numpart + ish - 1;
-           
-           if (fPyjets->K[0][in-1] == 1 || fPyjets->K[0][in-1] == 21 || fPyjets->K[0][in-1] == 11) {
-               jmin = in;
-               iGlu = in + 1;
-               iNew = in;
-           }
-           
-           kglu[i] = iGlu;
-//     
-// Shift stack
-//
-           for (Int_t j = jmax; j > jmin; j--)
-           {
-               for (Int_t k = 0; k < 5; k++) {
-                   fPyjets->K[k][j] =  fPyjets->K[k][j-ish];
-                   fPyjets->P[k][j] =  fPyjets->P[k][j-ish];
-                   fPyjets->V[k][j] =  fPyjets->V[k][j-ish];
-               }
-           } // end shifting
-           
+           iGlu = numpart;
+           if (iglu == 0) igMin = iGlu;
+           igMax = iGlu;
            numpart += ish;
            (fPyjets->N) += ish;
-               
+           
            if (ish == 1) {
-               fPyjets->P[0][iGlu] = p2[i][0];
-               fPyjets->P[1][iGlu] = p2[i][1];
-               fPyjets->P[2][iGlu] = p2[i][2];
-               fPyjets->P[3][iGlu] = p2[i][3];
-               fPyjets->P[4][iGlu] = p2[i][4];
+               fPyjets->P[0][iGlu] = p2[isys][0];
+               fPyjets->P[1][iGlu] = p2[isys][1];
+               fPyjets->P[2][iGlu] = p2[isys][2];
+               fPyjets->P[3][iGlu] = p2[isys][3];
+               fPyjets->P[4][iGlu] = p2[isys][4];
                
-               fPyjets->K[0][iGlu] = 2;
+               fPyjets->K[0][iGlu] = 1;
+               if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu] = 1;
                fPyjets->K[1][iGlu] = 21;       
-               fPyjets->K[2][iGlu] = fPyjets->K[2][iNew] + 1000;
+               fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
                fPyjets->K[3][iGlu] = -1;       
                fPyjets->K[4][iGlu] = -1;
+               
+               pg[0] += p2[isys][0];
+               pg[1] += p2[isys][1];
+               pg[2] += p2[isys][2];
+               pg[3] += p2[isys][3];
            } else {
                //
                // Split gluon in rest frame.
                //
-               Double_t bx   =  p2[i][0] / p2[i][3];
-               Double_t by   =  p2[i][1] / p2[i][3];
-               Double_t bz   =  p2[i][2] / p2[i][3];
-               Double_t pst  =  p2[i][4] / 2.;
+               Double_t bx   =  p2[isys][0] / p2[isys][3];
+               Double_t by   =  p2[isys][1] / p2[isys][3];
+               Double_t bz   =  p2[isys][2] / p2[isys][3];
+               Double_t pst  =  p2[isys][4] / 2.;
                //
                // Isotropic decay ????
                Double_t cost = 2. * gRandom->Rndm() - 1.;
@@ -1011,9 +966,9 @@ void  AliPythia::Quench()
                fPyjets->P[3][iGlu] = pst;
                fPyjets->P[4][iGlu] = 0.;
                
-               fPyjets->K[0][iGlu] = 2;
+               fPyjets->K[0][iGlu] = ;
                fPyjets->K[1][iGlu] = 21;       
-               fPyjets->K[2][iGlu] = fPyjets->K[2][iNew] + 1000;
+               fPyjets->K[2][iGlu] = fPyjets->K[2][in] + 1000;
                fPyjets->K[3][iGlu] = -1;       
                fPyjets->K[4][iGlu] = -1;
                
@@ -1023,9 +978,10 @@ void  AliPythia::Quench()
                fPyjets->P[3][iGlu+1] = pst;
                fPyjets->P[4][iGlu+1] = 0.;
                
-               fPyjets->K[0][iGlu+1] = 2;
+               fPyjets->K[0][iGlu+1] = 1;
+               if (iglu == nGluon[isys] - 1) fPyjets->K[0][iGlu+1] = 1;
                fPyjets->K[1][iGlu+1] = 21;     
-               fPyjets->K[2][iGlu+1] = fPyjets->K[2][iNew] + 1000;
+               fPyjets->K[2][iGlu+1] = fPyjets->K[2][in] + 1000;
                fPyjets->K[3][iGlu+1] = -1;     
                fPyjets->K[4][iGlu+1] = -1;
                SetMSTU(1,0);
@@ -1035,9 +991,35 @@ void  AliPythia::Quench()
                //
                Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
            }
-       }
+/*    
+           for (Int_t ig = iGlu; ig < iGlu+ish; ig++) {
+               Double_t px, py, pz;
+               px = fPyjets->P[0][ig]; 
+               py = fPyjets->P[1][ig]; 
+               pz = fPyjets->P[2][ig]; 
+               TVector3 v(px, py, pz);
+               v.RotateZ(-phiq[isys]);
+               v.RotateY(-thetaq[isys]);
+               Double_t pxs     = v.X(); Double_t pys = v.Y(); Double_t pzs  = v.Z();     
+               Double_t r       = AliPythiaRndm::GetPythiaRandom()->Rndm();
+               Double_t jtKick  = 0.3 * TMath::Sqrt(-TMath::Log(r));
+               if (ish == 2)   jtKick  = wjtKick[i] * TMath::Sqrt(-TMath::Log(r)) / TMath::Sqrt(2.);
+               Double_t phiKick = 2. * TMath::Pi() * AliPythiaRndm::GetPythiaRandom()->Rndm();
+               pxs += jtKick * TMath::Cos(phiKick);
+               pys += jtKick * TMath::Sin(phiKick);
+               TVector3 w(pxs, pys, pzs);
+               w.RotateY(thetaq[isys]);
+               w.RotateZ(phiq[isys]);
+               fPyjets->P[0][ig] = w.X(); 
+               fPyjets->P[1][ig] = w.Y(); 
+               fPyjets->P[2][ig] = w.Z(); 
+               fPyjets->P[2][ig] = w.Mag();
+           }
+*/
+       } // kGluon         
+       
        
-       // Check energy conservation
+    // Check energy conservation
        Double_t pxs = 0.;
        Double_t pys = 0.;
        Double_t pzs = 0.;      
@@ -1056,22 +1038,17 @@ void  AliPythia::Quench()
            TMath::Abs(pys) > 1.e-2 ||
            TMath::Abs(pzs) > 1.e-1) {
            printf("%e %e %e %e\n", pxs, pys, pzs, es);
-           this->Pylist(1);
-               Fatal("Quench()", "4-Momentum non-conservation");
+//             Fatal("Quench()", "4-Momentum non-conservation");
        }
-//     this->Pylist(1);
-    
-    } // end quenching loop
+       
+    } // end quenching loop (systems)
 // Clean-up
     for (Int_t i = 0; i < numpart; i++)
     {
-           imo =  fPyjets->K[2][i];
-           if (imo > 1000) fPyjets->K[2][i] = fPyjets->K[2][i] % 1000;
+       imo =  fPyjets->K[2][i];
+       if (imo > 1000) {
+           fPyjets->K[2][i] = fPyjets->K[2][i] % 1000;
+       }
     }
-    
-
-    
-//    delete[]    kNew;
-//    delete[]    pNew;
-    
+//     this->Pylist(1);
 } // end quench