Quenching: Problem with neg. inv. masses cured (z reduced for heavy quarks).
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 8 Apr 2004 16:54:30 +0000 (16:54 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 8 Apr 2004 16:54:30 +0000 (16:54 +0000)
PYTHIA6/AliPythia.cxx

index 0a8d044..85fd6d1 100644 (file)
@@ -655,9 +655,7 @@ void  AliPythia::Quench()
            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]);
        }
-       
        quenched[j] = (zInitial[j] > 0.01);
-       
     }
   
 //
@@ -669,261 +667,261 @@ void  AliPythia::Quench()
     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 k = 0; k < 4; k++)
-       {
-           p0[0][k] = 0.; p0[1][k] = 0.;
-           p1[0][k] = 0.; p1[1][k] = 0.;
-           p2[0][k] = 0.; p2[1][k] = 0.;
-       }
-       
-       Int_t nq[2] = {0, 0};
-       
-       for (Int_t i = 0; i < numpart; i++)
-       {
-           imo =  fPyjets->K[2][i];
-           kst =  fPyjets->K[0][i];
-           pdg =  fPyjets->K[1][i];
+       for (Int_t isys = 0; isys < 2; isys++) {
+//      Skip to next system if not quenched.
 
            
+           Double_t zHeavy = zInitial[isys];
            
+           if (!quenched[isys]) continue;
+//
+           while (1) {
+               icount = 0;
+               for (Int_t k = 0; k < 4; k++)
+               {
+                   p0[isys][k] = 0.;
+                   p1[isys][k] = 0.;
+                   p2[isys][k] = 0.;
+               }
+//      Loop over partons
+               for (Int_t i = 0; i < numpart; i++)
+               {
+                   imo =  fPyjets->K[2][i];
+                   kst =  fPyjets->K[0][i];
+                   pdg =  fPyjets->K[1][i];
+                   
+               
+               
 //      Quarks and gluons only
-           if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
+                   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 != 7 && imo != 8 && imo != 1007 && imo != 1008) continue;
-           
+                   if (imo > 8 && imo < 1000) imo = fPyjets->K[2][imo - 1];
+                   if (imo != isys + 7 && imo != 1000 + isys + 7) continue;
+                   
 //      Skip comment lines
-           if (kst != 1 && kst != 2) continue;
+                   if (kst != 1 && kst != 2) continue;
 //
 //      Parton kinematic
-           px    = fPyjets->P[0][i];
-           py    = fPyjets->P[1][i];
-           pz    = fPyjets->P[2][i];
-           e     = fPyjets->P[3][i];
-           m     = fPyjets->P[4][i];
-           pt    = TMath::Sqrt(px * px + py * py);
-           p     = TMath::Sqrt(px * px + py * py + pz * pz); 
-           phi   = TMath::Pi() + TMath::ATan2(-py, -px);
-           theta = TMath::ATan2(pt, pz);
-
+                   px    = fPyjets->P[0][i];
+                   py    = fPyjets->P[1][i];
+                   pz    = fPyjets->P[2][i];
+                   e     = fPyjets->P[3][i];
+                   m     = fPyjets->P[4][i];
+                   pt    = TMath::Sqrt(px * px + py * py);
+                   p     = TMath::Sqrt(px * px + py * py + pz * pz); 
+                   phi   = TMath::Pi() + TMath::ATan2(-py, -px);
+                   theta = TMath::ATan2(pt, pz);
+               
 //
 //      Save 4-momentum sum for balancing      
-           Int_t index = imo - 7;
-           if (index >=  1000) index -= 1000;
-
-           p0[index][0] += px;
-           p0[index][1] += py;
-           p0[index][2] += pz;
-           p0[index][3] += e;
-           
+                   Int_t index = imo - 7;
+                   if (index >=  1000) index -= 1000;
+                   
+                   p0[index][0] += px;
+                   p0[index][1] += py;
+                   p0[index][2] += pz;
+                   p0[index][3] += e;
+                   
 //      Don't quench radiated gluons
 //
-           if (imo == 1007 || imo == 1008) {
-               p1[index][0] += px;
-               p1[index][1] += py;
-               p1[index][2] += pz;
-               p1[index][3] += e;      
-               continue;
-           }
-           
+                   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;
+               
+                   klast[index] = i;
+                   
 //
 //      Fractional energy loss
-           Double_t z = zInitial[index];
-           if (!quenched[index]) continue;
-           //
-           //
-           //      Transform into frame in which initial parton is along z-axis
-           //
-           TVector3 v(px, py, pz);
-           v.RotateZ(-phiq[index]);
-           v.RotateY(-thetaq[index]);
-           Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pl  = v.Z();
-           Double_t jt  = TMath::Sqrt(pxs * pxs + pys * pys);
-           Double_t mt2 = jt * jt + m * m;
-           Double_t mt  = TMath::Sqrt(mt2);
-           Double_t zmin = 0.;
-           Double_t zmax = 1.;     
-           //
-           // z**2 mt**2 - pt**2 + pt'**2 > 0
-            // z**2 mt**2 + mt'**2 + m**2 - pt**2
-            // z**2 mt**2 + (1-z)**2 mt**2 + m**2 - pt**2
-            // mt**2(z**2 + 1 + z**2  - 2z) + m**2 - pt**2
-            // mt**2(2z**2 + 1 - 2z) + m**2 - pt**2 > 0
-            // mt**2(2z**2 + 1 - 2z) + 2 m**2 - mt**2 > 0
-            // mt**2(2z**2 - 2z) + 2 m**2 > 0
-            // z mt**2 (1 - z) -  m**2 < 0
-           // z**2 - z + 1/4 > 1/4 - m**2/mt**2
-            // (z-1/2)**2 > 1/4 - m**2/mt**2
-            // |z-1/2| > sqrt(1/4 - m**2/mt**2)
-            //
-            // m/mt < 1/2
-            // mt   > 2m
-           //
-           if (mt < 2. * m) {
-               printf("No phase space for quenching !: mt (%e) < 2 m (%e) \n", mt, m);
-               p1[index][0] += px;
-               p1[index][1] += py;
-               p1[index][2] += pz;
-               p1[index][3] += e;
-               continue;
-           } else {
-               zmin = 0.5 - TMath::Sqrt(0.25 - m * m / mt2);
-               if (z < zmin) {
-                   printf("No phase space for quenching ??: z (%e) < zmin (%e) \n", z, zmin);
-//                 z = zmin * 1.01;
-
-                   p1[index][0] += px;
-                   p1[index][1] += py;
-                   p1[index][2] += pz;
-                   p1[index][3] += e;
-                   continue;
-
-               }
-           }
-           //
-           // 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) {
+                   Double_t z = zInitial[index];
+                   if (m > 0.) z = zHeavy;
+                   
+                   //
+                   //
+                   //      Transform into frame in which initial parton is along z-axis
+                   //
+                   TVector3 v(px, py, pz);
+                   v.RotateZ(-phiq[index]);  v.RotateY(-thetaq[index]);
+                   Double_t pxs = v.X(); Double_t pys = v.Y(); Double_t pl  = v.Z();
+
+                   Double_t jt  = TMath::Sqrt(pxs * pxs + pys * pys);
+                   Double_t mt2 = jt * jt + m * m;
+                   Double_t zmax = 1.;     
+                   //
+                   // 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);
+                           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
                    
-                   p1[index][0] += px;
-                   p1[index][1] += py;
-                   p1[index][2] += pz;
-                   p1[index][3] += e;
-                   continue;
-               }
-           } // massive particles
-           
-           //
-           // Change light-cone kinematics rel. to initial parton
-           //  
-           Double_t eppzOld = e + pl;
-           Double_t empzOld = e - pl;
-           
-           Double_t eppzNew = (1. - z) * eppzOld;
-           Double_t empzNew = empzOld - mt2 * z / eppzOld;
-           Double_t eNew0    = 0.5 * (eppzNew + empzNew);
-           Double_t pzNew0   = 0.5 * (eppzNew - empzNew);
-           
-           Double_t ptNew;
-           //
-           // 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.;
+                   //
+                   // Change light-cone kinematics rel. to initial parton
+                   //  
+                   Double_t eppzOld = e + pl;
+                   Double_t empzOld = e - pl;
+                   
+                   Double_t eppzNew = (1. - z) * eppzOld;
+                   Double_t empzNew = empzOld - mt2 * z / eppzOld;
+                   Double_t eNew    = 0.5 * (eppzNew + empzNew);
+                   Double_t plNew   = 0.5 * (eppzNew - empzNew);
+                   
+                   Double_t jtNew;
+                   //
+                   // 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;
+                   } else {
+                       jtNew    = TMath::Sqrt(mt2New - m * m);
+                   }
+               
            
-           if (m * m > mt2New) {
+                   //
+                   //     Calculate new px, py
+                   //
+                   Double_t pxNew   = jtNew / jt * pxs;
+                   Double_t pyNew   = jtNew / jt * pys;        
+                   
+//                 Double_t dpx = pxs - pxNew;
+//                 Double_t dpy = pys - pyNew;
+//                 Double_t dpz = pl  - plNew;
+//                 Double_t de  = e   - eNew;
+//                 Double_t dmass2 = de * de  - dpx * dpx - dpy * dpy - dpz * dpz;
+//                 printf("New mass (1) %e %e %e %e %e %e %e \n", dmass2, jt, jtNew, pl, plNew, e, eNew);
+//                 printf("New mass (2) %e %e \n", pxNew, pyNew);
+                   //
+                   //      Rotate back
+                   //  
+                   TVector3 w(pxNew, pyNew, plNew);
+                   w.RotateY(thetaq[index]); w.RotateZ(phiq[index]);
+                   pxNew = w.X(); pyNew = w.Y(); plNew = w.Z();
+               
+                   p1[index][0] += pxNew;
+                   p1[index][1] += pyNew;
+                   p1[index][2] += plNew;
+                   p1[index][3] += eNew;       
+                   //
+                   // Updated 4-momentum vectors
+                   //
+                   pNew[icount][0]  = pxNew;
+                   pNew[icount][1]  = pyNew;
+                   pNew[icount][2]  = plNew;
+                   pNew[icount][3]  = eNew;
+                   kNew[icount]     = i;
+                   icount++;
+               } // parton loop
                //
-               // This should not happen 
+               // Check if there was phase-space for quenching
                //
-               Fatal("Quench()", "This should never happen %e %e %e!", m, eppzNew, empzNew);
-               ptNew = 0;
-           } else {
-               ptNew    = TMath::Sqrt(mt2New - m * m);
-           }
 
+               if (icount == 0) quenched[isys] = kFALSE;
+               if (!quenched[isys]) break;
+               
+               for (Int_t j = 0; j < 4; j++) 
+               {
+                   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]);
+                   if (p2[isys][4] < -0.01) {
+                       printf("Negative mass squared ! Let's try to fix this by decreasing z\n");
+//                     this->Pylist(1);
+                       
+                   } 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) {
+                   printf("No success ! \n");
+                   icount = 0;
+                   quenched[isys] = kFALSE;
+                   break;
+               }
+           } // iteration on z
            
-           //
-           //     Calculate new px, py
-           //
-           Double_t pxNew0   = ptNew / jt * pxs;
-           Double_t pyNew0   = ptNew / jt * pys;       
-/*
-           Double_t dpx = pxs - pxNew0;
-           Double_t dpy = pys - pyNew0;
-           Double_t dpz = pl  - pzNew0;
-           Double_t de  = e   - eNew0;
-           Double_t dmass2 = de * de  - dpx * dpx - dpy * dpy - dpz * dpz;
-*/
-           //
-           //      Rotate back
-           //  
-           TVector3 w(pxNew0, pyNew0, pzNew0);
-           w.RotateY(thetaq[index]);
-           w.RotateZ(phiq[index]);
-           pxNew0 = w.X(); pyNew0 = w.Y(); pzNew0 = w.Z();
-           
-           
-           p1[index][0] += pxNew0;
-           p1[index][1] += pyNew0;
-           p1[index][2] += pzNew0;
-           p1[index][3] += eNew0;      
-           //
-           // Update event record
-           //
-           fPyjets->P[0][i] = pxNew0;
-           fPyjets->P[1][i] = pyNew0;
-           fPyjets->P[2][i] = pzNew0;
-           fPyjets->P[3][i] = eNew0;
-           nq[index]++;
-           
-       }
-       
-       //
-       // Gluons
-       // 
-       
-       for (Int_t k = 0; k < 2; k++) 
-       {
-           //
-           // Check if there was phase-space for quenching
-           //
-           if (nq[k] == 0) quenched[k] = kFALSE;
-           
-           if (!quenched[k]) continue;
-           
-           for (Int_t j = 0; j < 4; j++) 
-           {
-               p2[k][j] = p0[k][j] - p1[k][j];
+//         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] );
+               fPyjets->P[0][kNew[k]] = pNew[k][0];
+               fPyjets->P[1][kNew[k]] = pNew[k][1];
+               fPyjets->P[2][kNew[k]] = pNew[k][2];
+               fPyjets->P[3][kNew[k]] = pNew[k][3];
            }
-           p2[k][4] = p2[k][3] * p2[k][3] - p2[k][0] * p2[k][0] - p2[k][1] * p2[k][1] - p2[k][2] * p2[k][2];
-
-           if (p2[k][4] > 0.) {
-               p2[k][4] = TMath::Sqrt(p2[k][4]);
-           } else {
-               printf("Warning negative mass squared in system %d %f ! \n", k, zInitial[k]);
-               printf("Kinematics %10.3e %10.3e %10.3e %10.3e %10.3e \n", p2[k][0], p2[k][1], p2[k][2], p2[k][3], p2[k][4]);
-               if (p2[k][4] < -0.1) Fatal("Boost", "Negative mass squared !");
-               p2[k][4] = 0.;
-           }
-           //
-           // 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]);
-           */
-       }
-       
+       } // System         
        //
        // Add the gluons
        //
-       
        Int_t ish = 0;    
        for (Int_t i = 0; i < 2; i++) {
            Int_t jmin, jmax, iGlu, iNew;
@@ -971,7 +969,7 @@ void  AliPythia::Quench()
            
            numpart += ish;
            (fPyjets->N) += ish;
-           
+               
            if (ish == 1) {
                fPyjets->P[0][iGlu] = p2[i][0];
                fPyjets->P[1][iGlu] = p2[i][1];
@@ -1037,8 +1035,8 @@ void  AliPythia::Quench()
                //
                Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
            }
-       } // end adding gluons
-       //
+       }
+       
        // Check energy conservation
        Double_t pxs = 0.;
        Double_t pys = 0.;
@@ -1059,15 +1057,21 @@ void  AliPythia::Quench()
            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");
        }
-
-    } // end quenchin loop
-    // Clean-up
+//     this->Pylist(1);
+    
+    } // end quenching loop
+// Clean-up
     for (Int_t i = 0; i < numpart; i++)
     {
-       imo =  fPyjets->K[2][i];
-       if (imo > 1000) 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;
+    
 } // end quench