]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PYTHIA6/AliPythia.cxx
Including of Lambda(1520) to the generation. Small fraction of Lambda0 (0.2), Sigma0...
[u/mrichter/AliRoot.git] / PYTHIA6 / AliPythia.cxx
index f83a094bcfd5444968cbd72e2496bb29c5682408..22b27e7c0450b7385a95f9e1a16286a80410035d 100644 (file)
@@ -28,11 +28,15 @@ ClassImp(AliPythia)
 # define pycell pycell_
 # define pyshow pyshow_
 # define pyrobo pyrobo_
+# define pyquen pyquen_
+# define pyevnw pyevnw_
 # define type_of_call
 #else
 # define pyclus PYCLUS
 # define pycell PYCELL
 # define pyrobo PYROBO
+# define pyquen PYQUEN
+# define pyevnw PYEVNW
 # define type_of_call _stdcall
 #endif
 
@@ -40,6 +44,8 @@ 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 &);
+extern "C" void type_of_call pyquen(Double_t &, Int_t &, Double_t &);
+extern "C" void type_of_call pyevnw(){;}
 
 //_____________________________________________________________________________
 
@@ -65,9 +71,17 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
     fProcess = process;
     fEcms = energy;
     fStrucFunc = strucfunc;
-//  don't decay p0
-    SetMDCY(Pycomp(111),1,0);
-//  select structure function 
+//...Switch off decay of pi0, K0S, Lambda, Sigma+-, Xi0-, Omega-.
+    SetMDCY(Pycomp(111) ,1,0);
+    SetMDCY(Pycomp(310) ,1,0);
+    SetMDCY(Pycomp(3122),1,0);
+    SetMDCY(Pycomp(3112),1,0);
+    SetMDCY(Pycomp(3212),1,0);
+    SetMDCY(Pycomp(3222),1,0);
+    SetMDCY(Pycomp(3312),1,0);
+    SetMDCY(Pycomp(3322),1,0);
+    SetMDCY(Pycomp(3334),1,0);
+    //  select structure function 
     SetMSTP(52,2);
     SetMSTP(51,strucfunc);
 //
@@ -81,9 +95,71 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
 //  select charm production
     switch (process) 
     {
+    case kPyOldUEQ2ordered:  //Old underlying events with Q2 ordered QCD processes
+//        Multiple interactions on.
+       SetMSTP(81,1);
+// Double Gaussian matter distribution.
+       SetMSTP(82,4);
+       SetPARP(83,0.5);
+       SetPARP(84,0.4);
+//  pT0.
+       SetPARP(82,2.0);
+//  Reference energy for pT0 and energy rescaling pace.
+       SetPARP(89,1800);
+       SetPARP(90,0.25);
+//  String drawing almost completely minimizes string length.
+       SetPARP(85,0.9);
+       SetPARP(86,0.95);
+// ISR and FSR activity.
+       SetPARP(67,4);
+       SetPARP(71,4);
+// Lambda_FSR scale.
+       SetPARJ(81,0.29);
+       break;
+    case kPyOldUEQ2ordered2:   
+// Old underlying events with Q2 ordered QCD processes
+// Multiple interactions on.
+       SetMSTP(81,1);
+// Double Gaussian matter distribution.
+       SetMSTP(82,4);
+       SetPARP(83,0.5);
+       SetPARP(84,0.4);
+// pT0.
+       SetPARP(82,2.0);
+// Reference energy for pT0 and energy rescaling pace.
+       SetPARP(89,1800);
+       SetPARP(90,0.16);  // here is the difference with  kPyOldUEQ2ordered
+// String drawing almost completely minimizes string length.
+       SetPARP(85,0.9);
+       SetPARP(86,0.95);
+// ISR and FSR activity.
+       SetPARP(67,4);
+       SetPARP(71,4);
+// Lambda_FSR scale.
+       SetPARJ(81,0.29);       
+       break;
+    case kPyOldPopcorn:  
+// Old production mechanism: Old Popcorn
+       SetMSEL(1);
+       SetMSTJ(12,3); 
+// (D=2) Like MSTJ(12)=2 but added prod ofthe 1er rank baryon
+       SetMSTP(88,2); 
+// (D=1)see can be used to form  baryons (BARYON JUNCTION)
+       SetMSTJ(1,1);  
+       SetMSTP(51,kCTEQ5L);// CTEQ 5L        ! CTEQ5L pdf
+       SetMSTP(81,1);      // Multiple Interactions ON
+       SetMSTP(82,4);      // Double Gaussian Model         
+       SetPARP(82,1.8);    // [GeV]    PT_min at Ref. energy
+       SetPARP(89,1000.);  // [GeV]   Ref. energy
+       SetPARP(90,0.16);   // 2*epsilon (exponent in power law)
+       SetPARP(83,0.5);    // Core density in proton matter dist. (def.value)
+       SetPARP(84,0.5);    // Core radius
+       SetPARP(85,0.33);   // Regulates gluon prod. mechanism
+       SetPARP(86,0.66);   // Regulates gluon prod. mechanism
+       SetPARP(67,1);      // Regulate gluon prod. mechanism
+       break;
     case kPyCharm:
        SetMSEL(4);
-//
 //  heavy quark masses
 
        SetPMAS(4,1,1.2);
@@ -174,7 +250,7 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
 // ATLAS Tuning
 //
        
-       SetMSTP(51,7);             // CTEQ5L pdf
+       SetMSTP(51,kCTEQ5L);       // CTEQ5L pdf
        SetMSTP(81,1);             // Multiple Interactions ON
        SetMSTP(82,4);             // Double Gaussian Model
 
@@ -192,7 +268,17 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
 //  QCD Jets
 //
        SetMSEL(1);
-       break;
+ // Pythia Tune A (CDF)
+ //
+       SetPARP(67,4.);            // Regulates Initial State Radiation
+       SetMSTP(82,4);             // Double Gaussian Model
+       SetPARP(82,2.0);           // [GeV]    PT_min at Ref. energy
+       SetPARP(84,0.4);           // Core radius
+       SetPARP(85,0.90) ;         // Regulates gluon prod. mechanism
+       SetPARP(86,0.95);          // Regulates gluon prod. mechanism
+       SetPARP(89,1800.);         // [GeV]   Ref. energy
+       SetPARP(90,0.25);          // 2*epsilon (exponent in power law)
+       break;
     case kPyDirectGamma:
        SetMSEL(10);
        break;
@@ -538,6 +624,55 @@ void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfun
       SetPMAS(5,1,4.75);
 
       break;
+    case kPyW:
+
+      //Inclusive production of W+/-
+      SetMSEL(0);
+      //f fbar -> W+ 
+      SetMSUB(2,1);
+      //       //f fbar -> g W+
+      //       SetMSUB(16,1);
+      //       //f fbar -> gamma W+
+      //       SetMSUB(20,1);
+      //       //f g -> f W+  
+      //       SetMSUB(31,1);
+      //       //f gamma -> f W+
+      //       SetMSUB(36,1);
+      
+      // Initial/final parton shower on (Pythia default)
+      // With parton showers on we are generating "W inclusive process"
+      SetMSTP(61,1); //Initial QCD & QED showers on
+      SetMSTP(71,1); //Final QCD & QED showers on
+      
+      break;  
+
+    case kPyZ:
+
+      //Inclusive production of Z
+      SetMSEL(0);
+      //f fbar -> Z/gamma
+      SetMSUB(1,1);
+      
+      //       // f fbar -> g Z/gamma
+      //       SetMSUB(15,1);
+      //       // f fbar -> gamma Z/gamma
+      //       SetMSUB(19,1);
+      //       // f g -> f Z/gamma
+      //       SetMSUB(30,1);
+      //       // f gamma -> f Z/gamma
+      //       SetMSUB(35,1);
+      
+      //only Z included, not gamma
+      SetMSTP(43,2);
+      
+      // Initial/final parton shower on (Pythia default)
+      // With parton showers on we are generating "Z inclusive process"
+      SetMSTP(61,1); //Initial QCD & QED showers on
+      SetMSTP(71,1); //Final QCD & QED showers on
+      
+      break;  
+
     }
 //
 //  Initialize PYTHIA
@@ -703,10 +838,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
 //
@@ -756,9 +892,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.
@@ -769,7 +905,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 ++;
@@ -781,37 +917,41 @@ 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
 
            
-           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]);
+//         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]));
 
 
@@ -878,7 +1018,7 @@ void  AliPythia::Quench()
                    
 //
 //      Fractional energy loss
-                   Double_t z = zInitial[index];
+                   Double_t z = zquench[index];
                    
                    
 //      Don't fully quench radiated gluons
@@ -993,7 +1133,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");
@@ -1191,3 +1331,26 @@ void  AliPythia::Quench()
 //     this->Pylist(1);
 } // end quench
 
+
+void AliPythia::Pyquen(Double_t a, Int_t ibf, Double_t b)
+{
+    // Igor Lokthine's quenching routine
+    pyquen(a, ibf, b);
+}
+
+void AliPythia::Pyevnw()
+{
+    // New multiple interaction scenario
+    pyevnw();
+}
+
+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];
+
+}
+
+