Simple jet quenching algorythm added.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 24 Sep 2003 12:34:38 +0000 (12:34 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 24 Sep 2003 12:34:38 +0000 (12:34 +0000)
PYTHIA6/AliPythia.cxx
PYTHIA6/AliPythia.h

index 5613d528c7b47c0850316f23c198efcfaad9acca..798e76a9e8230c73aece941d82d2dfb1cbdaaa0d 100644 (file)
@@ -23,15 +23,20 @@ ClassImp(AliPythia)
 #ifndef WIN32
 # define pyclus pyclus_
 # define pycell pycell_
 #ifndef WIN32
 # define pyclus pyclus_
 # define pycell pycell_
+# define pyshow pyshow_
+# define pyrobo pyrobo_
 # define type_of_call
 #else
 # define pyclus PYCLUS
 # define pycell PYCELL
 # define type_of_call
 #else
 # define pyclus PYCLUS
 # define pycell PYCELL
+# define pyrobo PYROBO
 # define type_of_call _stdcall
 #endif
 
 extern "C" void type_of_call pyclus(Int_t & );
 extern "C" void type_of_call pycell(Int_t & );
 # define type_of_call _stdcall
 #endif
 
 extern "C" void type_of_call pyclus(Int_t & );
 extern "C" void type_of_call pycell(Int_t & );
+extern "C" void type_of_call pyshow(Int_t &, Int_t &, Double_t &);
+extern "C" void type_of_call pyrobo(Int_t &, Int_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
 
 //_____________________________________________________________________________
 
 
 //_____________________________________________________________________________
 
@@ -497,3 +502,288 @@ void  AliPythia::Pycell(Int_t& njet)
     pycell(njet);
 }
 
     pycell(njet);
 }
 
+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);
+}
+
+void AliPythia::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez)
+{
+    pyrobo(imi, ima, the, phi, bex, bey, bez);
+}
+
+
+
+void  AliPythia::Quench()
+{
+//
+//
+//  Simple Jet Quenching routine:
+//  =============================
+//  The jet formed by all final state partons radiated by the parton created 
+//  in the hard collisions is quenched by a factor z using:
+//  (E + p_z)new = (1-z) (E + p_z)old
+//
+//   The lost momentum is first balanced by one gluon with virtuality > 0.   
+//   Subsequently the gluon splits to yield two gluons with E = p.
+//
+    Float_t p0[2][5];
+    Float_t p1[2][5];
+    Float_t p2[2][5];
+    Int_t   klast[2] = {-1, -1};
+    Int_t   kglu[2];
+    for (Int_t i = 0; i < 4; i++)
+    {
+       p0[0][i] = 0.;
+       p0[1][i] = 0.;
+       p1[0][i] = 0.;
+       p1[1][i] = 0.;
+       p2[0][i] = 0.;
+       p2[1][i] = 0.;
+    }
+
+    Int_t numpart   = fPyjets->N;
+
+    for (Int_t i = 0; i < numpart; i++)
+    {
+       Int_t imo =  fPyjets->K[2][i];
+       Int_t kst =  fPyjets->K[0][i];
+       Int_t pdg =  fPyjets->K[1][i];
+       
+//      Quarks and gluons only
+       if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
+
+//      Particles from hard scattering only
+
+
+       Float_t px = fPyjets->P[0][i];
+       Float_t py = fPyjets->P[1][i];
+       Float_t pz = fPyjets->P[2][i];
+       Float_t e  = fPyjets->P[3][i];
+       Float_t m  = fPyjets->P[4][i];
+       Float_t pt = TMath::Sqrt(px * px + py * py);
+//      Skip comment lines
+       if (kst != 1 && kst != 2) continue;
+
+       Float_t mt = TMath::Sqrt(px * px + py * py + m * m);
+
+       //
+       // Some cuts to be in a save kinematic region
+       //
+       if (imo != 7 && imo != 8) continue;
+       Int_t index = imo - 7;
+       klast[index] = i;
+       
+       p0[index][0] += px;
+       p0[index][1] += py;
+       p0[index][2] += pz;
+       p0[index][3] += e;      
+       
+       //
+       // Fix z
+       //
+
+       Float_t z       = 0.2;
+       Float_t eppzOld = e + pz;
+       Float_t empzOld = e - pz;
+
+
+       //
+       // Kinematics of the original parton
+       // 
+
+       Float_t eppzNew = (1. - z) * eppzOld;
+       Float_t empzNew = empzOld - mt * mt * z / eppzOld;
+       Float_t eNew0    = 0.5 * (eppzNew + empzNew);
+       Float_t pzNew0   = 0.5 * (eppzNew - empzNew);
+       //
+        // Skip if pt too small
+       // 
+       if (m * m > eppzNew * empzNew) continue;
+       Float_t ptNew    = TMath::Sqrt(eppzNew * empzNew - m * m);
+       Float_t pxNew0   = ptNew / pt * px;
+       Float_t pyNew0   = ptNew / pt * py;     
+
+       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;
+       
+    }
+
+    //
+    // Gluons
+    // 
+    
+    for (Int_t k = 0; k < 2; k++) 
+    {
+       for (Int_t j = 0; j < 4; j++) 
+       {
+           p2[k][j] = p0[k][j] - p1[k][j];
+       }
+       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.)
+       {
+
+           //
+           // Bring gluon back to mass shell via momentum scaling 
+           // (momentum will not be conserved, but energy)
+           //
+           // not used anymore
+/*
+           Float_t psq   = p2[k][0] * p2[k][0] + p2[k][1] * p2[k][1] + p2[k][2] * p2[k][2];
+           Float_t fact  = TMath::Sqrt(1. + p2[k][4] / psq);
+           p2[k][0] *= fact;
+           p2[k][1] *= fact;
+           p2[k][2] *= fact;
+           p2[k][3]  = TMath::Sqrt(psq) * fact;
+           p2[k][4]  = 0.;
+*/
+       }
+    }
+
+    if (p2[0][4] > 0.) {
+       p2[0][4] = TMath::Sqrt(p2[0][4]);
+    } else {
+       printf("Warning negative mass squared ! \n");
+    }
+    
+    if (p2[1][4] > 0.) {
+       p2[1][4] = TMath::Sqrt(p2[1][4]);
+    } else {
+       printf("Warning negative mass squared ! \n");
+    }
+    
+    //
+    // Add the gluons
+    //
+    
+    
+    for (Int_t i = 0; i < 2; i++) {
+       Int_t ish, jmin, jmax, iGlu, iNew;
+       Int_t in = klast[i];
+       ish = 0;
+       
+       if (in == -1) continue;
+       if (i == 1 && klast[1] > klast[0]) in += ish;
+       
+       jmin = in - 1;
+       ish  = 1;
+
+       if (p2[i][4] > 0) ish = 2;
+       
+       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;
+       
+       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
+       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->K[0][iGlu] = 2;
+           fPyjets->K[1][iGlu] = 21;   
+           fPyjets->K[2][iGlu] = fPyjets->K[2][iNew];
+           fPyjets->K[3][iGlu] = -1;   
+           fPyjets->K[4][iGlu] = -1;
+       } 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];
+           
+           Float_t pst  = p2[i][4] / 2.;
+           
+           Float_t cost = 2. * gRandom->Rndm() - 1.;
+           Float_t sint = TMath::Sqrt(1. - cost * cost);
+           Float_t phi = 2. * TMath::Pi() * gRandom->Rndm();
+           
+           Float_t pz1 =   pst * cost;
+           Float_t pz2 =  -pst * cost;
+           Float_t pt1 =   pst * sint;
+           Float_t pt2 =  -pst * sint;
+           Float_t px1 =   pt1 * TMath::Cos(phi);
+           Float_t py1 =   pt1 * TMath::Sin(phi);          
+           Float_t px2 =   pt2 * TMath::Cos(phi);
+           Float_t py2 =   pt2 * TMath::Sin(phi);          
+           
+           fPyjets->P[0][iGlu] = px1;
+           fPyjets->P[1][iGlu] = py1;
+           fPyjets->P[2][iGlu] = pz1;
+           fPyjets->P[3][iGlu] = pst;
+           fPyjets->P[4][iGlu] = 0.;
+           
+           fPyjets->K[0][iGlu] = 2;
+           fPyjets->K[1][iGlu] = 21;   
+           fPyjets->K[2][iGlu] = fPyjets->K[2][iNew];
+           fPyjets->K[3][iGlu] = -1;   
+           fPyjets->K[4][iGlu] = -1;
+
+           fPyjets->P[0][iGlu+1] = px2;
+           fPyjets->P[1][iGlu+1] = py2;
+           fPyjets->P[2][iGlu+1] = pz2;
+           fPyjets->P[3][iGlu+1] = pst;
+           fPyjets->P[4][iGlu+1] = 0.;
+           
+           fPyjets->K[0][iGlu+1] = 2;
+           fPyjets->K[1][iGlu+1] = 21; 
+           fPyjets->K[2][iGlu+1] = fPyjets->K[2][iNew];
+           fPyjets->K[3][iGlu+1] = -1; 
+           fPyjets->K[4][iGlu+1] = -1;
+           SetMSTU(1,0);
+           SetMSTU(2,0);
+
+           //
+           // Boost back
+           //
+           Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
+
+       }
+    } // end adding gluons
+} // end quench
+
+
index b0d7375b017e3cffa1c63914d0106e573c5a0c74..c91fa9872ef2c545618ea34012865901cd4e8cb6 100644 (file)
@@ -51,9 +51,12 @@ class AliPythia : public TPythia6, public AliRndm
     virtual void SetDecayTable();
     virtual void Pycell(Int_t& nclus);
     virtual void Pyclus(Int_t& nclus);
     virtual void SetDecayTable();
     virtual void Pycell(Int_t& nclus);
     virtual void Pyclus(Int_t& nclus);
+    virtual void Pyshow(Int_t ip1, Int_t ip2, Double_t qmax);
+    virtual void Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez);
+    
     // return instance of the singleton
     static  AliPythia* Instance();
     // return instance of the singleton
     static  AliPythia* Instance();
-
+    virtual void Quench();
  protected:
     Process_t     fProcess;           // Process type
     Float_t       fEcms;              // Centre of mass energy
  protected:
     Process_t     fProcess;           // Process type
     Float_t       fEcms;              // Centre of mass energy