]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PYTHIA6/AliGenPythiaPlus.cxx
make ptmin setteable
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythiaPlus.cxx
index 5579ebd84f138816b35bc03238936a36beacb28e..0f91b97fcce8f5933e52f5a29aea32c3347d784c 100644 (file)
@@ -124,8 +124,9 @@ AliGenPythiaPlus::AliGenPythiaPlus():
     fPHOSEta(0.13),
     fEMCALMinPhi(79.),
     fEMCALMaxPhi(191.),
-    fEMCALEta(0.71)
-
+    fEMCALEta(0.71),
+    fItune(-1), 
+    fInfo(1) 
 {
 // Default Constructor
   fEnergyCMS = 5500.;
@@ -211,7 +212,9 @@ AliGenPythiaPlus::AliGenPythiaPlus(AliPythiaBase* pythia)
      fPHOSEta(0.13),
      fEMCALMinPhi(79.),
      fEMCALMaxPhi(191.),
-     fEMCALEta(0.71)
+     fEMCALEta(0.71),
+     fItune(-1),
+     fInfo(1) 
 {
 // default charm production at 5. 5 TeV
 // semimuonic decay
@@ -345,7 +348,7 @@ void AliGenPythiaPlus::Init()
        fRL = 0x0;
     }
  //
-    fPythia->ProcInit(fProcess, fEnergyCMS, fStrucFunc);
+    fPythia->ProcInit(fProcess, fEnergyCMS, fStrucFunc, fItune);
     //  Forward Paramters to the AliPythia object
     fDecayer->SetForceDecay(fForceDecay);    
 // Switch off Heavy Flavors on request  
@@ -371,6 +374,7 @@ void AliGenPythiaPlus::Init()
     case kPyCharmpPbMNR:
     case kPyCharmppMNR:
     case kPyCharmppMNRwmi:
+    case kPyCharmPWHG:
        fParentSelect[0] =   411;
        fParentSelect[1] =   421;
        fParentSelect[2] =   431;
@@ -398,11 +402,17 @@ void AliGenPythiaPlus::Init()
        fParentSelect[0] =   431;
        fFlavorSelect    =   4; 
        break;
+    case kPyLambdacppMNR:
+       fParentSelect[0] =  4122;
+       fFlavorSelect    =   4; 
+       break;      
     case kPyBeauty:
+    case kPyBeautyJets:
     case kPyBeautyPbPbMNR:
     case kPyBeautypPbMNR:
     case kPyBeautyppMNR:
     case kPyBeautyppMNRwmi:
+    case kPyBeautyPWHG:
        fParentSelect[0]=  511;
        fParentSelect[1]=  521;
        fParentSelect[2]=  531;
@@ -426,17 +436,24 @@ void AliGenPythiaPlus::Init()
     case kPyJpsi:
        fParentSelect[0] = 443;
        break;
+    case kPyMbAtlasTuneMC09:
     case kPyMbDefault:
     case kPyMb:
     case kPyMbWithDirectPhoton:
     case kPyMbNonDiffr:
     case kPyMbMSEL1:
     case kPyJets:
+    case kPyJetsPWHG:
     case kPyDirectGamma:
     case kPyLhwgMb:    
        break;
+    case kPyWPWHG:
     case kPyW:
     case kPyZ:
+    case kPyZgamma:
+    case kPyMBRSingleDiffraction:
+    case kPyMBRDoubleDiffraction:
+    case kPyMBRCentralDiffraction:
         break;
     }
 //
@@ -491,9 +508,9 @@ void AliGenPythiaPlus::Generate()
     
     fDecayer->ForceDecay();
 
-    Float_t polar[3]   =   {0,0,0};
-    Float_t origin[3]  =   {0,0,0};
-    Float_t p[4];
+    Double_t polar[3]   =   {0,0,0};
+    Double_t origin[3]  =   {0,0,0};
+    Double_t p[4];
 //  converts from mm/c to s
     const Float_t kconv=0.001/2.999792458e8;
 //
@@ -529,8 +546,8 @@ void AliGenPythiaPlus::Generate()
            }
            fNpartons = fPythia->GetNumberOfParticles();
        } else {
-           printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
-           fRL->GetEvent(AliRunLoader::GetRunLoader()->GetEventNumber());
+           printf("Loading Event %d\n",AliRunLoader::Instance()->GetEventNumber());
+           fRL->GetEvent(AliRunLoader::Instance()->GetEventNumber());
            fPythia->SetNumberOfParticles(0);
            fPythia->LoadEvent(fRL->Stack(), 0 , 1);
            fPythia->EditEventList(21);
@@ -591,8 +608,13 @@ void AliGenPythiaPlus::Generate()
            fProcess != kPyMbMSEL1     &&
            fProcess != kPyW && 
            fProcess != kPyZ &&
+      fProcess != kPyZgamma &&
            fProcess != kPyCharmppMNRwmi && 
-           fProcess != kPyBeautyppMNRwmi) {
+           fProcess != kPyBeautyppMNRwmi &&
+      fProcess != kPyWPWHG &&
+           fProcess != kPyJetsPWHG &&
+           fProcess != kPyCharmPWHG &&
+     fProcess != kPyBeautyPWHG) {
            
            for (i = 0; i < np; i++) {
                TParticle* iparticle = (TParticle *) fParticles.At(i);
@@ -742,7 +764,7 @@ void AliGenPythiaPlus::Generate()
                    origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
                    origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
                    
-                   Float_t tof   = kconv*iparticle->T();
+                   Float_t tof   = fTime + kconv*iparticle->T();
                    Int_t ipa = (fPythia->Version() == 6) ? (iparticle->GetFirstMother() - 1) :(iparticle->GetFirstMother()) ;
                    Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
  
@@ -783,7 +805,7 @@ void AliGenPythiaPlus::Generate()
          }
            if (jev >= fNpart || fNpart == -1) {
                fKineBias=Float_t(fNpart)/Float_t(fTrials);
-               fPythia->GetXandQ(fQ, fX1, fX2);
+               if (fInfo) fPythia->GetXandQ(fX1, fX2, fQ);
                fTrialsRun += fTrials;
                fNev++;
                MakeHeader();
@@ -815,7 +837,7 @@ Int_t  AliGenPythiaPlus::GenerateMB()
     
     Int_t* pParent = new Int_t[np];
     for (i=0; i< np; i++) pParent[i] = -1;
-    if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
+    if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyJetsPWHG || fProcess == kPyCharmPWHG || fProcess == kPyBeautyPWHG ) {
        TParticle* jet1 = (TParticle *) fParticles.At(6);
        TParticle* jet2 = (TParticle *) fParticles.At(7);
        if (!CheckTrigger(jet1, jet2)) {
@@ -825,7 +847,7 @@ Int_t  AliGenPythiaPlus::GenerateMB()
     }
 
     // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
-    if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
+    if ((fProcess == kPyJets || fProcess == kPyJetsPWHG) && (fFragPhotonInCalo || fPi0InCalo) ) {
 
       Bool_t ok = kFALSE;
 
@@ -850,13 +872,15 @@ Int_t  AliGenPythiaPlus::GenerateMB()
            }
        }
       }
-      if(!ok)
-       return 0;
+      if(!ok){
+         delete [] pParent;
+         return 0;
+      }
     }
     
     
      // Select events with a photon  pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
-    if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
+    if ((fProcess == kPyJets || fProcess == kPyJetsPWHG || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
 
       Bool_t okd = kFALSE;
 
@@ -885,8 +909,10 @@ Int_t  AliGenPythiaPlus::GenerateMB()
       if(!okd && iphcand != -1) // execute rotation in phi 
           RotatePhi(iphcand,okd);
       
-      if(!okd)
-       return 0;
+      if(!okd) {
+         delete[] pParent;
+         return 0;
+      }
     }
     
     if (fTriggerParticle) {
@@ -931,7 +957,7 @@ Int_t  AliGenPythiaPlus::GenerateMB()
          if(mi<0) continue;
          mother = (TParticle*)fParticles.At(mi);
          mpdg=TMath::Abs(mother->GetPdgCode());
-         mpdgUpperFamily=(mpdg>1000 ? mpdg-1000 : mpdg-100); // keep e from c from b
+         mpdgUpperFamily=(mpdg>1000 ? mpdg+1000 : mpdg+100); // keep e from c from b
          if ( ParentSelected(mpdg) || 
              (fFlavorSelect==5 && ParentSelected(mpdgUpperFamily))) {
            if (KinematicSelection(partCheck,1)) {
@@ -952,8 +978,11 @@ Int_t  AliGenPythiaPlus::GenerateMB()
     }
 
     //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
-    if ( (fProcess == kPyW ||
+    if ( (
+    fProcess == kPyWPWHG ||
+    fProcess == kPyW ||
          fProcess == kPyZ ||
+    fProcess == kPyZgamma ||
          fProcess == kPyMbDefault ||
          fProcess == kPyMb ||
          fProcess == kPyMbWithDirectPhoton ||
@@ -974,7 +1003,7 @@ Int_t  AliGenPythiaPlus::GenerateMB()
        Int_t km = iparticle->GetFirstMother();
        if ((ks == 1  && kf!=0 && KinematicSelection(iparticle, 0)) ||
            (ks != 1) ||
-           (fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
+           ((fProcess == kPyJets || fProcess == kPyJetsPWHG) && ks == 21 && km == 0 && i>1)) {
            nc++;
            if (ks == 1) trackIt = 1;
 
@@ -994,7 +1023,7 @@ Int_t  AliGenPythiaPlus::GenerateMB()
            origin[1] = fVertex[1]+iparticle->Vy()/10; // [cm]
            origin[2] = fVertex[2]+iparticle->Vz()/10; // [cm]
            
-           Float_t tof = fEventTime + kconv * iparticle->T();
+           Float_t tof = fTime + fEventTime + kconv * iparticle->T();
 
            PushTrack(fTrackIt*trackIt, iparent, kf, 
                      p[0], p[1], p[2], p[3], 
@@ -1007,17 +1036,18 @@ Int_t  AliGenPythiaPlus::GenerateMB()
            //
            // Special Treatment to store color-flow
            //
+           /*
            if (ks == 3 || ks == 13 || ks == 14) {
                TParticle* particle = 0;
                if (fStack) {
                    particle = fStack->Particle(nt);
                } else {
-                   particle = gAlice->Stack()->Particle(nt);
+                   particle = AliRunLoader::Instance()->Stack()->Particle(nt);
                }
-//             particle->SetFirstDaughter(fPythia->GetK(2, i));
-//             particle->SetLastDaughter(fPythia->GetK(3, i));         
+               particle->SetFirstDaughter(fPythia->GetK(2, i));
+               particle->SetLastDaughter(fPythia->GetK(3, i));         
            }
-           
+           */  
            KeepTrack(nt);
            pParent[i] = nt;
            SetHighWaterMark(nt);
@@ -1092,14 +1122,14 @@ void AliGenPythiaPlus::MakeHeader()
 //
 // Event Vertex 
     fHeader->SetPrimaryVertex(fVertex);
-    
+    fHeader->SetInteractionTime(fTime+fEventTime);    
 //
 // Number of primaries
     fHeader->SetNProduced(fNprimaries);
 //
 // Jets that have triggered
 
-    if (fProcess == kPyJets)
+    if (fProcess == kPyJets || fProcess == kPyJetsPWHG)
     {
        Int_t ntrig, njet;
        Float_t jets[4][10];
@@ -1131,8 +1161,9 @@ void AliGenPythiaPlus::MakeHeader()
 // Store quenching parameters
 //
     if (fQuench){
-       Double_t z[4];
-       Double_t xp, yp;
+        Double_t z[4] = {0.};
+       Double_t xp = 0.;
+       Double_t yp = 0.;
        if (fQuench == 1) {
            // Pythia::Quench()
            fPythia->GetQuenchingParameters(xp, yp, z);
@@ -1160,7 +1191,7 @@ void AliGenPythiaPlus::MakeHeader()
     fHeader = 0x0;
 }
 
-Bool_t AliGenPythiaPlus::CheckTrigger(TParticle* jet1, TParticle* jet2)
+Bool_t AliGenPythiaPlus::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
 {
 // Check the kinematic trigger condition
 //
@@ -1175,7 +1206,7 @@ Bool_t AliGenPythiaPlus::CheckTrigger(TParticle* jet1, TParticle* jet2)
     pdg[1] = jet2->GetPdgCode();    
     Bool_t   triggered = kFALSE;
 
-    if (fProcess == kPyJets) {
+    if (fProcess == kPyJets || fProcess == kPyJetsPWHG) {
        Int_t njets = 0;
        Int_t ntrig = 0;
        Float_t jets[4][10];
@@ -1329,7 +1360,7 @@ void AliGenPythiaPlus::GetSubEventTime()
 
 
 
-Bool_t AliGenPythiaPlus::IsInEMCAL(Float_t phi, Float_t eta)
+Bool_t AliGenPythiaPlus::IsInEMCAL(Float_t phi, Float_t eta) const
 {
   // Is particle in EMCAL acceptance? 
   // phi in degrees, etamin=-etamax
@@ -1340,7 +1371,7 @@ Bool_t AliGenPythiaPlus::IsInEMCAL(Float_t phi, Float_t eta)
     return kFALSE;
 }
 
-Bool_t AliGenPythiaPlus::IsInPHOS(Float_t phi, Float_t eta)
+Bool_t AliGenPythiaPlus::IsInPHOS(Float_t phi, Float_t eta) const
 {
   // Is particle in PHOS acceptance? 
   // Acceptance slightly larger considered.