X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=PYTHIA6%2FAliGenPythia.cxx;h=d60a161704caa863979ce88d269350cd39902272;hp=7b65e67151399bd613355ff27cf271a3ec21a650;hb=38db0ee672885e3cc0687e0a2519da99846dd843;hpb=f913ec4fd8eb28e70b23988bfbb1c3628465697c diff --git a/PYTHIA6/AliGenPythia.cxx b/PYTHIA6/AliGenPythia.cxx index 7b65e671513..d60a161704c 100644 --- a/PYTHIA6/AliGenPythia.cxx +++ b/PYTHIA6/AliGenPythia.cxx @@ -54,6 +54,10 @@ AliGenPythia::AliGenPythia() fHeader = 0; fReadFromFile = 0; fEventTime = 0.; + fInteractionRate = 0.; + fTimeWindow = 0.; + fEventsTime = 0; + fCurSubEvent = 0; fDecayer = new AliDecayerPythia(); SetEventListRange(); SetJetPhiRange(); @@ -64,7 +68,10 @@ AliGenPythia::AliGenPythia() SetPtKick(); SetQuench(); SetHadronisation(); + SetTriggerParticle(); fSetNuclei = kFALSE; + fNewMIS = kFALSE; + fHFoff = kFALSE; if (!AliPythiaRndm::GetPythiaRandom()) AliPythiaRndm::SetPythiaRandom(GetRandom()); } @@ -81,6 +88,10 @@ AliGenPythia::AliGenPythia(Int_t npart) fXsection = 0.; fReadFromFile = 0; fEventTime = 0.; + fInteractionRate = 0.; + fTimeWindow = 0.; + fEventsTime = 0; + fCurSubEvent = 0; SetProcess(); SetStrucFunc(); SetForceDecay(); @@ -105,6 +116,7 @@ AliGenPythia::AliGenPythia(Int_t npart) 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 @@ -115,6 +127,8 @@ AliGenPythia::AliGenPythia(Int_t npart) // Pycel SetPycellParameters(); fSetNuclei = kFALSE; + fNewMIS = kFALSE; + fHFoff = kFALSE; } AliGenPythia::AliGenPythia(const AliGenPythia & Pythia) @@ -127,6 +141,57 @@ AliGenPythia::AliGenPythia(const AliGenPythia & Pythia) 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, @@ -202,19 +267,28 @@ 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); // Parent and Children Selection switch (fProcess) { + case kPyOldUEQ2ordered: + case kPyOldUEQ2ordered2: + case kPyOldPopcorn: + break; case kPyCharm: case kPyCharmUnforced: case kPyCharmPbPbMNR: - case kPyCharmppMNR: case kPyCharmpPbMNR: + case kPyCharmppMNR: + case kPyCharmppMNRwmi: fParentSelect[0] = 411; fParentSelect[1] = 421; fParentSelect[2] = 431; @@ -233,10 +307,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; @@ -262,10 +343,12 @@ void AliGenPythia::Init() break; case kPyMb: case kPyMbNonDiffr: + case kPyMbMSEL1: case kPyJets: case kPyDirectGamma: break; case kPyW: + case kPyZ: break; } // @@ -274,17 +357,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; @@ -303,11 +386,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() @@ -347,7 +429,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()); @@ -403,7 +489,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); @@ -592,8 +680,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); @@ -629,6 +716,7 @@ 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) { @@ -636,7 +724,54 @@ Int_t AliGenPythia::GenerateMB() 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; iAt(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 && yQAt(i); @@ -694,7 +829,6 @@ Int_t AliGenPythia::GenerateMB() if (pParent) delete[] pParent; - printf("\n I've put %i particles on the stack \n",nc); return 1; } @@ -811,13 +945,25 @@ void AliGenPythia::MakeHeader() ((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) { @@ -866,6 +1012,36 @@ Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2) } 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; jAt(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) { @@ -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 - 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