]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PYTHIA6/AliGenPythia.cxx
Changed inheritance order (RenderElement base comes first).
[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx
index c6405e2b906a188d0d1308edb2f3acb903c2c041..ddbbe4fd0c89489f86b065b2356ddad1048495cb 100644 (file)
@@ -68,7 +68,14 @@ AliGenPythia::AliGenPythia()
   SetPtKick();
   SetQuench();
   SetHadronisation();  
+  SetTriggerParticle();
+  SetNuclei(0,0);
   fSetNuclei = kFALSE;
+  fNewMIS    = kFALSE;
+  fHFoff     = kFALSE;
+  fGinit     = 1;
+  fGfinal    = 1;
+  
   if (!AliPythiaRndm::GetPythiaRandom()) 
     AliPythiaRndm::SetPythiaRandom(GetRandom());
 }
@@ -113,6 +120,8 @@ AliGenPythia::AliGenPythia(Int_t npart)
     SetQuench();
     SetHadronisation();
     SetPtKick();
+    SetTriggerParticle();
+    SetNuclei(0,0);
     // Options determining what to keep in the stack (Heavy flavour generation)
     fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
     fFeedDownOpt = kTRUE;             // allow feed down from higher family
@@ -123,7 +132,11 @@ AliGenPythia::AliGenPythia(Int_t npart)
     // Pycel
     SetPycellParameters();
     fSetNuclei = kFALSE;
-}
+    fNewMIS    = kFALSE;
+    fHFoff     = kFALSE;
+    fGinit     = 1;
+    fGfinal    = 1;
+ }
 
 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
     :AliGenMC(Pythia)
@@ -261,8 +274,12 @@ void AliGenPythia::Init()
     } else {
        fRL = 0x0;
     }
-    
-    
+// Switch off Heavy Flavors on request  
+    if (fHFoff) {
+       fPythia->SetMSTP(58, 3);
+       fPythia->SetMSTJ(45, 3);        
+       for (Int_t i = 156; i <= 160; i++) fPythia->SetMDME(i, 1, 0);
+    }
  //
     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
 
@@ -276,8 +293,9 @@ void AliGenPythia::Init()
     case kPyCharm:
     case kPyCharmUnforced:
     case kPyCharmPbPbMNR:
-    case kPyCharmppMNR:
     case kPyCharmpPbMNR:
+    case kPyCharmppMNR:
+    case kPyCharmppMNRwmi:
        fParentSelect[0] =   411;
        fParentSelect[1] =   421;
        fParentSelect[2] =   431;
@@ -296,10 +314,17 @@ void AliGenPythia::Init()
        fParentSelect[0] =   411;
        fFlavorSelect    =   4; 
        break;
+    case kPyDPlusStrangePbPbMNR:
+    case kPyDPlusStrangepPbMNR:
+    case kPyDPlusStrangeppMNR:
+       fParentSelect[0] =   431;
+       fFlavorSelect    =   4; 
+       break;
     case kPyBeauty:
     case kPyBeautyPbPbMNR:
     case kPyBeautypPbMNR:
     case kPyBeautyppMNR:
+    case kPyBeautyppMNRwmi:
        fParentSelect[0]=  511;
        fParentSelect[1]=  521;
        fParentSelect[2]=  531;
@@ -325,10 +350,12 @@ void AliGenPythia::Init()
        break;
     case kPyMb:
     case kPyMbNonDiffr:
+    case kPyMbMSEL1:
     case kPyJets:
     case kPyDirectGamma:
        break;
     case kPyW:
+    case kPyZ:
         break;
     }
 //
@@ -337,17 +364,17 @@ void AliGenPythia::Init()
 //
 //  Configure detector (EMCAL like)
 //
-       fPythia->SetPARU(51, fPycellEtaMax);
-       fPythia->SetMSTU(51, fPycellNEta);
-       fPythia->SetMSTU(52, fPycellNPhi);
+    fPythia->SetPARU(51, fPycellEtaMax);
+    fPythia->SetMSTU(51, fPycellNEta);
+    fPythia->SetMSTU(52, fPycellNPhi);
 //
 //  Configure Jet Finder
 //  
-       fPythia->SetPARU(58,  fPycellThreshold);
-       fPythia->SetPARU(52,  fPycellEtSeed);
-       fPythia->SetPARU(53,  fPycellMinEtJet);
-       fPythia->SetPARU(54,  fPycellMaxRadius);
-       fPythia->SetMSTU(54,  2);
+    fPythia->SetPARU(58,  fPycellThreshold);
+    fPythia->SetPARU(52,  fPycellEtSeed);
+    fPythia->SetPARU(53,  fPycellMinEtJet);
+    fPythia->SetPARU(54,  fPycellMaxRadius);
+    fPythia->SetMSTU(54,  2);
 //
 //  This counts the total number of calls to Pyevnt() per run.
     fTrialsRun = 0;
@@ -366,11 +393,10 @@ void AliGenPythia::Init()
        fDyBoost = 0;
        Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
     }
-
+    
     if (fQuench) {
        fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
     }
-    
 }
 
 void AliGenPythia::Generate()
@@ -410,7 +436,11 @@ void AliGenPythia::Generate()
 // Either produce new event or read partons from file
 //     
        if (!fReadFromFile) {
-           fPythia->Pyevnt();
+           if (!fNewMIS) {
+               fPythia->Pyevnt();
+           } else {
+               fPythia->Pyevnw();
+           }
            fNpartons = fPythia->GetN();
        } else {
            printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());
@@ -466,7 +496,9 @@ void AliGenPythia::Generate()
        if (fProcess != kPyMb && fProcess != kPyJets && 
            fProcess != kPyDirectGamma &&
            fProcess != kPyMbNonDiffr  &&
-           fProcess != kPyW) {
+           fProcess != kPyMbMSEL1     &&
+           fProcess != kPyW && fProcess != kPyZ &&
+           fProcess != kPyCharmppMNRwmi && fProcess != kPyBeautyppMNRwmi) {
            
            for (i = 0; i < np; i++) {
                TParticle* iparticle = (TParticle *) fParticles->At(i);
@@ -634,9 +666,9 @@ void AliGenPythia::Generate()
        
        GetSubEventTime();
 
-       if (pParent)   delete[] pParent;
-       if (pSelected) delete[] pSelected;
-       if (trackIt)   delete[] trackIt;
+       delete[] pParent;
+       delete[] pSelected;
+       delete[] trackIt;
 
        if (nc > 0) {
          switch (fCountMode) {
@@ -655,8 +687,7 @@ void AliGenPythia::Generate()
          }
            if (jev >= fNpart || fNpart == -1) {
                fKineBias=Float_t(fNpart)/Float_t(fTrials);
-               printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
-
+               
                fQ  += fPythia->GetVINT(51);
                fX1 += fPythia->GetVINT(41);
                fX2 += fPythia->GetVINT(42);
@@ -692,18 +723,65 @@ Int_t  AliGenPythia::GenerateMB()
     Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
 
 
+
     Int_t* pParent = new Int_t[np];
     for (i=0; i< np; i++) pParent[i] = -1;
     if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
        TParticle* jet1 = (TParticle *) fParticles->At(6);
        TParticle* jet2 = (TParticle *) fParticles->At(7);
-       if (!CheckTrigger(jet1, jet2)) return 0;
+       if (!CheckTrigger(jet1, jet2)) {
+         delete [] pParent;
+         return 0;
+       }
+    }
+
+    if (fTriggerParticle) {
+       Bool_t triggered = kFALSE;
+       for (i = 0; i < np; i++) {
+           TParticle *  iparticle = (TParticle *) fParticles->At(i);
+           kf = CheckPDGCode(iparticle->GetPdgCode());
+           if (kf != fTriggerParticle) continue;
+           if (iparticle->Pt() == 0.) continue;
+           if (TMath::Abs(iparticle->Eta()) > fTriggerEta) continue;
+           triggered = kTRUE;
+           break;
+       }
+       if (!triggered) {
+         delete [] pParent;
+         return 0;
+       }
+    }
+       
+
+    // Check if there is a ccbar or bbbar pair with at least one of the two
+    // in fYMin < y < fYMax
+    if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi) {
+      TParticle *hvq;
+      Bool_t  theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
+      Float_t yQ;  
+      Int_t   pdgQ;
+      for(i=0; i<np; i++) {
+       hvq = (TParticle*)fParticles->At(i);
+       pdgQ = hvq->GetPdgCode();  
+       if(TMath::Abs(pdgQ) != fFlavorSelect) continue; 
+       if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
+       yQ = 0.5*TMath::Log((hvq->Energy()+hvq->Pz()+1.e-13)/
+                           (hvq->Energy()-hvq->Pz()+1.e-13));
+       if(yQ>fYMin && yQ<fYMax) inYcut=kTRUE;
+      }
+      if (!theQ || !theQbar || !inYcut) {
+       delete[] pParent;
+       return 0;
+      }
     }
-    
 
-    //Introducing child cuts in case kPyW
-    if ( (fProcess == kPyW)  && (fCutOnChild == 1) ) {
-       if ( !CheckKinematicsOnChild() ) return 0;
+    //Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
+    if ( (fProcess == kPyW || fProcess == kPyZ || fProcess == kPyMb || fProcess == kPyMbNonDiffr)  
+        && (fCutOnChild == 1) ) {
+      if ( !CheckKinematicsOnChild() ) {
+       delete[] pParent;
+       return 0;
+      }
     }
   
 
@@ -762,9 +840,8 @@ Int_t  AliGenPythia::GenerateMB()
        } // select particle
     } // particle loop 
 
-    if (pParent) delete[] pParent;
+    delete[] pParent;
     
-    printf("\n I've put %i particles on the stack \n",nc);
     return 1;
 }
 
@@ -784,7 +861,7 @@ void AliGenPythia::FinishRun()
     printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
 }
 
-void AliGenPythia::AdjustWeights()
+void AliGenPythia::AdjustWeights() const
 {
 // Adjust the weights after generation of all events
 //
@@ -810,6 +887,9 @@ void AliGenPythia::SetNuclei(Int_t a1, Int_t a2)
 
 void AliGenPythia::MakeHeader()
 {
+//
+// Make header for the simulated event
+// 
   if (gAlice) {
     if (gAlice->GetEvNumber()>=fDebugEventFirst &&
        gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(2);
@@ -881,7 +961,9 @@ void AliGenPythia::MakeHeader()
            ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
            ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
        }
-    
+//
+// Store pt^hard 
+    ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));
 //
 //  Pass header
 //
@@ -923,7 +1005,7 @@ Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
 //
        GetJets(njets, ntrig, jets);
        
-       if (ntrig) triggered = kTRUE;
+       if (ntrig || fEtMinJet == 0.) triggered = kTRUE;
 //
     } else {
        Int_t ij = 0;
@@ -948,9 +1030,11 @@ Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
 }
 
 
-//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
-Bool_t AliGenPythia::CheckKinematicsOnChild(){
 
+Bool_t AliGenPythia::CheckKinematicsOnChild(){
+//
+//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
+//
     Bool_t checking = kFALSE;
     Int_t j, kcode, ks, km;
     Int_t nPartAcc = 0; //number of particles in the acceptance range
@@ -958,7 +1042,7 @@ Bool_t AliGenPythia::CheckKinematicsOnChild(){
     if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
     Int_t npart = fParticles->GetEntriesFast();
     
-    for (j = 0; j<npart; j++) {  
+    for (j = 0; j<npart; j++) {
        TParticle *  jparticle = (TParticle *) fParticles->At(j);
        kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
        ks = jparticle->GetStatusCode();
@@ -967,14 +1051,13 @@ Bool_t AliGenPythia::CheckKinematicsOnChild(){
        if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
            nPartAcc++;
        }
+       if( numberOfAcceptedParticles <= nPartAcc){
+         checking = kTRUE;
+         break;
+       }
     }
-  
-    if( numberOfAcceptedParticles <= nPartAcc){
-       checking = kTRUE;
-    }
-    
+
     return checking;
-    
 }
 
          
@@ -1003,12 +1086,12 @@ void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
     
     
     for (Int_t part = 0; part < npart; part++) {
-       TParticle *MPart = stack->Particle(part);
+       TParticle *mPart = stack->Particle(part);
        
-       Int_t kf     =  MPart->GetPdgCode();
-       Int_t ks     =  MPart->GetStatusCode();
-       Int_t idf    =  MPart->GetFirstDaughter();
-       Int_t idl    =  MPart->GetLastDaughter();
+       Int_t kf     =  mPart->GetPdgCode();
+       Int_t ks     =  mPart->GetStatusCode();
+       Int_t idf    =  mPart->GetFirstDaughter();
+       Int_t idl    =  mPart->GetLastDaughter();
        
        if (reHadr) {
            if (ks == 11 || ks == 12) {
@@ -1018,11 +1101,11 @@ void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
            }
        }
        
-       Float_t px = MPart->Px();
-       Float_t py = MPart->Py();
-       Float_t pz = MPart->Pz();
-       Float_t e  = MPart->Energy();
-       Float_t m  = MPart->GetCalcMass();
+       Float_t px = mPart->Px();
+       Float_t py = mPart->Py();
+       Float_t pz = mPart->Pz();
+       Float_t e  = mPart->Energy();
+       Float_t m  = mPart->GetCalcMass();
        
        
        (fPythia->GetPyjets())->P[0][part+n0] = px;
@@ -1035,7 +1118,7 @@ void  AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
        (fPythia->GetPyjets())->K[0][part+n0] = ks;
        (fPythia->GetPyjets())->K[3][part+n0] = idf + 1;
        (fPythia->GetPyjets())->K[4][part+n0] = idl + 1;
-       (fPythia->GetPyjets())->K[2][part+n0] = MPart->GetFirstMother() + 1;
+       (fPythia->GetPyjets())->K[2][part+n0] = mPart->GetFirstMother() + 1;
     }
 }