Trigger on particle and antiparticle separately.

[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythia.cxx

diff --git a/PYTHIA6/AliGenPythia.cxx b/PYTHIA6/AliGenPythia.cxx
index 7b65e67151399bd613355ff27cf271a3ec21a650..d60a161704caa863979ce88d269350cd39902272 100644 (file)
--- a/PYTHIA6/AliGenPythia.cxx
+++ b/PYTHIA6/AliGenPythia.cxx
@@ -54,6 +54,10 @@ AliGenPythia::AliGenPythia()
   fHeader = 0;
   fReadFromFile = 0;
   fEventTime = 0.;
   fHeader = 0;
   fReadFromFile = 0;
   fEventTime = 0.;

+  fInteractionRate = 0.;
+  fTimeWindow = 0.;
+  fEventsTime = 0;
+  fCurSubEvent = 0;

   fDecayer   = new AliDecayerPythia();
   SetEventListRange();
   SetJetPhiRange();
   fDecayer   = new AliDecayerPythia();
   SetEventListRange();
   SetJetPhiRange();


@@ -64,7 +68,10 @@ AliGenPythia::AliGenPythia()
   SetPtKick();
   SetQuench();
   SetHadronisation();  
   SetPtKick();
   SetQuench();
   SetHadronisation();  

+  SetTriggerParticle();

   fSetNuclei = kFALSE;
   fSetNuclei = kFALSE;

+  fNewMIS    = kFALSE;
+  fHFoff     = kFALSE;

   if (!AliPythiaRndm::GetPythiaRandom()) 
     AliPythiaRndm::SetPythiaRandom(GetRandom());
 }
   if (!AliPythiaRndm::GetPythiaRandom()) 
     AliPythiaRndm::SetPythiaRandom(GetRandom());
 }


@@ -81,6 +88,10 @@ AliGenPythia::AliGenPythia(Int_t npart)
     fXsection  = 0.;
     fReadFromFile = 0;
     fEventTime = 0.;
     fXsection  = 0.;
     fReadFromFile = 0;
     fEventTime = 0.;

+    fInteractionRate = 0.;
+    fTimeWindow = 0.;
+    fEventsTime = 0;
+    fCurSubEvent = 0;

     SetProcess();
     SetStrucFunc();
     SetForceDecay();
     SetProcess();
     SetStrucFunc();
     SetForceDecay();


@@ -105,6 +116,7 @@ AliGenPythia::AliGenPythia(Int_t npart)
     SetQuench();
     SetHadronisation();
     SetPtKick();
     SetQuench();
     SetHadronisation();
     SetPtKick();

+    SetTriggerParticle();

     // 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
     // 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


@@ -115,6 +127,8 @@ AliGenPythia::AliGenPythia(Int_t npart)
     // Pycel
     SetPycellParameters();
     fSetNuclei = kFALSE;
     // Pycel
     SetPycellParameters();
     fSetNuclei = kFALSE;

+    fNewMIS    = kFALSE;
+    fHFoff     = kFALSE;

 }
 
 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
 }
 
 AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)


@@ -127,6 +141,57 @@ AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
 AliGenPythia::~AliGenPythia()
 {
 // Destructor
 AliGenPythia::~AliGenPythia()
 {
 // Destructor

+  if(fEventsTime) delete fEventsTime;
+}
+
+void AliGenPythia::SetInteractionRate(Float_t rate,Float_t timewindow)
+{
+// Generate pileup using user specified rate
+    fInteractionRate = rate;
+    fTimeWindow = timewindow;
+    GeneratePileup();
+}
+
+void AliGenPythia::GeneratePileup()
+{
+// Generate sub events time for pileup
+    fEventsTime = 0;
+    if(fInteractionRate == 0.) {
+      Warning("GeneratePileup","Zero interaction specified. Skipping pileup generation.\n");
+      return;
+    }
+
+    Int_t npart = NumberParticles();
+    if(npart < 0) {
+      Warning("GeneratePileup","Negative number of particles. Skipping pileup generation.\n");
+      return;
+    }
+
+    if(fEventsTime) delete fEventsTime;
+    fEventsTime = new TArrayF(npart);
+    TArrayF &array = *fEventsTime;
+    for(Int_t ipart = 0; ipart < npart; ipart++)
+      array[ipart] = 0.;
+
+    Float_t eventtime = 0.;
+    while(1)
+      {
+       eventtime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
+       if(eventtime > fTimeWindow) break;
+       array.Set(array.GetSize()+1);
+       array[array.GetSize()-1] = eventtime;
+      }
+
+    eventtime = 0.;
+    while(1)
+      {
+       eventtime -= (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
+       if(TMath::Abs(eventtime) > fTimeWindow) break;
+       array.Set(array.GetSize()+1);
+       array[array.GetSize()-1] = eventtime;
+      }
+
+    SetNumberParticles(fEventsTime->GetSize());

 }
 
 void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,
 }
 
 void AliGenPythia::SetPycellParameters(Float_t etamax, Int_t neta, Int_t nphi,


@@ -202,19 +267,28 @@ void AliGenPythia::Init()
     } else {
        fRL = 0x0;
     }
     } 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);
 
 //  Parent and Children Selection
     switch (fProcess) 
     {
  //
     fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
 
 //  Parent and Children Selection
     switch (fProcess) 
     {

+    case kPyOldUEQ2ordered:
+    case kPyOldUEQ2ordered2:
+    case kPyOldPopcorn:
+      break;

     case kPyCharm:
     case kPyCharmUnforced:
     case kPyCharmPbPbMNR:
     case kPyCharm:
     case kPyCharmUnforced:
     case kPyCharmPbPbMNR:

-    case kPyCharmppMNR:

     case kPyCharmpPbMNR:
     case kPyCharmpPbMNR:

+    case kPyCharmppMNR:
+    case kPyCharmppMNRwmi:

        fParentSelect[0] =   411;
        fParentSelect[1] =   421;
        fParentSelect[2] =   431;
        fParentSelect[0] =   411;
        fParentSelect[1] =   421;
        fParentSelect[2] =   431;


@@ -233,10 +307,17 @@ void AliGenPythia::Init()
        fParentSelect[0] =   411;
        fFlavorSelect    =   4; 
        break;
        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 kPyBeauty:
     case kPyBeautyPbPbMNR:
     case kPyBeautypPbMNR:
     case kPyBeautyppMNR:

+    case kPyBeautyppMNRwmi:

        fParentSelect[0]=  511;
        fParentSelect[1]=  521;
        fParentSelect[2]=  531;
        fParentSelect[0]=  511;
        fParentSelect[1]=  521;
        fParentSelect[2]=  531;


@@ -262,10 +343,12 @@ void AliGenPythia::Init()
        break;
     case kPyMb:
     case kPyMbNonDiffr:
        break;
     case kPyMb:
     case kPyMbNonDiffr:

+    case kPyMbMSEL1:

     case kPyJets:
     case kPyDirectGamma:
        break;
     case kPyW:
     case kPyJets:
     case kPyDirectGamma:
        break;
     case kPyW:

+    case kPyZ:

         break;
     }
 //
         break;
     }
 //


@@ -274,17 +357,17 @@ void AliGenPythia::Init()
 //
 //  Configure detector (EMCAL like)
 //
 //
 //  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
 //  
 //
 //  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;
 //
 //  This counts the total number of calls to Pyevnt() per run.
     fTrialsRun = 0;


@@ -303,11 +386,10 @@ void AliGenPythia::Init()
        fDyBoost = 0;
        Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n");
     }
        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);
     }
     if (fQuench) {
        fPythia->InitQuenching(0., 0.1, 0.6e6, 0);
     }

-    

 }
 
 void AliGenPythia::Generate()
 }
 
 void AliGenPythia::Generate()


@@ -347,7 +429,11 @@ void AliGenPythia::Generate()
 // Either produce new event or read partons from file
 //     
        if (!fReadFromFile) {
 // 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());
            fNpartons = fPythia->GetN();
        } else {
            printf("Loading Event %d\n",AliRunLoader::GetRunLoader()->GetEventNumber());


@@ -403,7 +489,9 @@ void AliGenPythia::Generate()
        if (fProcess != kPyMb && fProcess != kPyJets && 
            fProcess != kPyDirectGamma &&
            fProcess != kPyMbNonDiffr  &&
        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);
            
            for (i = 0; i < np; i++) {
                TParticle* iparticle = (TParticle *) fParticles->At(i);


@@ -592,8 +680,7 @@ void AliGenPythia::Generate()
          }
            if (jev >= fNpart || fNpart == -1) {
                fKineBias=Float_t(fNpart)/Float_t(fTrials);
          }
            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);
                fQ  += fPythia->GetVINT(51);
                fX1 += fPythia->GetVINT(41);
                fX2 += fPythia->GetVINT(42);


@@ -629,6 +716,7 @@ Int_t  AliGenPythia::GenerateMB()
     Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
 
 
     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) {
     Int_t* pParent = new Int_t[np];
     for (i=0; i< np; i++) pParent[i] = -1;
     if (fProcess == kPyJets || fProcess == kPyDirectGamma) {


@@ -636,7 +724,54 @@ Int_t  AliGenPythia::GenerateMB()
        TParticle* jet2 = (TParticle *) fParticles->At(7);
        if (!CheckTrigger(jet1, jet2)) return 0;
     }
        TParticle* jet2 = (TParticle *) fParticles->At(7);
        if (!CheckTrigger(jet1, jet2)) 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) 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) {
+       if (pParent) delete[] pParent;
+       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() ) {
+       if (pParent)   delete[] pParent;
+       return 0;
+      }
+    }
+  
+

     for (i = 0; i < np; i++) {
        Int_t trackIt = 0;
        TParticle *  iparticle = (TParticle *) fParticles->At(i);
     for (i = 0; i < np; i++) {
        Int_t trackIt = 0;
        TParticle *  iparticle = (TParticle *) fParticles->At(i);


@@ -694,7 +829,6 @@ Int_t  AliGenPythia::GenerateMB()
 
     if (pParent) delete[] pParent;
     
 
     if (pParent) delete[] pParent;
     

-    printf("\n I've put %i particles on the stack \n",nc);

     return 1;
 }
 
     return 1;
 }
 


@@ -811,13 +945,25 @@ void AliGenPythia::MakeHeader()
            ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
            ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
        }
            ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp);
            ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z);
        }

-    

 //
 //

-//  Pass header to RunLoader
+// Store pt^hard 
+    ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47));

 //
 //

-    AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(fHeader);   
+//  Pass header
+//
+    AddHeader(fHeader);

 }
 }

-       
+
+void AliGenPythia::AddHeader(AliGenEventHeader* header)
+{
+    // Add header to container or runloader
+    if (fContainer) {
+       fContainer->AddHeader(header);
+    } else {
+       AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);   
+    }
+}
+

 
 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
 {
 
 Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
 {


@@ -866,6 +1012,36 @@ Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
     }
     return triggered;
 }
     }
     return triggered;
 }

+
+
+//Checking Kinematics on Child (status code 1, particle code ?, kin cuts
+Bool_t AliGenPythia::CheckKinematicsOnChild(){
+
+    Bool_t checking = kFALSE;
+    Int_t j, kcode, ks, km;
+    Int_t nPartAcc = 0; //number of particles in the acceptance range
+    Int_t numberOfAcceptedParticles = 1;
+    if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
+    Int_t npart = fParticles->GetEntriesFast();
+    
+    for (j = 0; j<npart; j++) {
+       TParticle *  jparticle = (TParticle *) fParticles->At(j);
+       kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
+       ks = jparticle->GetStatusCode();
+       km = jparticle->GetFirstMother(); 
+       
+       if( (ks == 1)  &&  (kcode == fPdgCodeParticleforAcceptanceCut)  &&  (KinematicSelection(jparticle,1)) ){
+           nPartAcc++;
+       }
+       if( numberOfAcceptedParticles <= nPartAcc){
+         checking = kTRUE;
+         break;
+       }
+    }
+
+    return checking;
+}
+

          
 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
 {
          
 AliGenPythia& AliGenPythia::operator=(const  AliGenPythia& rhs)
 {


@@ -1018,7 +1194,13 @@ void  AliGenPythia::GetJets(Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10])
 void AliGenPythia::GetSubEventTime()
 {
   // Calculates time of the next subevent
 void AliGenPythia::GetSubEventTime()
 {
   // Calculates time of the next subevent

-  fEventTime += (AliPythiaRndm::GetPythiaRandom())->Exp(1./fInteractionRate);
+  fEventTime = 0.;
+  if (fEventsTime) {
+    TArrayF &array = *fEventsTime;
+    fEventTime = array[fCurSubEvent++];
+  }
+  //  printf(" Event time: %d %f %p",fCurSubEvent,fEventTime,fEventsTime);
+  return;

 }
 
 #ifdef never
 }
 
 #ifdef never