X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=PYTHIA6%2FAliGenPythia.cxx;h=8566ef5885fa5969cf43921de4bf4b7af8d8fc94;hb=4ec650ac1817219804a6fa716bc842e93a263e22;hp=83cbdd02893e4d7fd19634174fbe85aa574b1e72;hpb=1058d9df13456512344aae22e9b810e032811eab;p=u%2Fmrichter%2FAliRoot.git diff --git a/PYTHIA6/AliGenPythia.cxx b/PYTHIA6/AliGenPythia.cxx index 83cbdd02893..8566ef5885f 100644 --- a/PYTHIA6/AliGenPythia.cxx +++ b/PYTHIA6/AliGenPythia.cxx @@ -35,6 +35,7 @@ #include "AliConst.h" #include "AliDecayerPythia.h" #include "AliGenPythia.h" +#include "AliFastGlauber.h" #include "AliHeader.h" #include "AliGenPythiaEventHeader.h" #include "AliPythia.h" @@ -51,6 +52,7 @@ ClassImp(AliGenPythia) AliGenPythia::AliGenPythia(): AliGenMC(), fProcess(kPyCharm), + fItune(-1), fStrucFunc(kCTEQ5L), fKineBias(0.), fTrials(0), @@ -74,9 +76,13 @@ AliGenPythia::AliGenPythia(): fGinit(1), fGfinal(1), fHadronisation(1), + fPatchOmegaDalitz(0), fNpartons(0), fReadFromFile(0), fQuench(0), + fQhat(0.), + fLength(0.), + fImpact(0.), fPtKick(1.), fFullEvent(kTRUE), fDecayer(new AliDecayerPythia()), @@ -111,18 +117,21 @@ AliGenPythia::AliGenPythia(): fTriggerEta(0.9), fTriggerMultiplicity(0), fTriggerMultiplicityEta(0), + fTriggerMultiplicityPtMin(0), fCountMode(kCountAll), fHeader(0), fRL(0), - fFileName(0), + fkFileName(0), fFragPhotonInCalo(kFALSE), fPi0InCalo(kFALSE) , fPhotonInCalo(kFALSE), + fEleInEMCAL(kFALSE), fCheckEMCAL(kFALSE), fCheckPHOS(kFALSE), fCheckPHOSeta(kFALSE), fFragPhotonOrPi0MinPt(0), fPhotonMinPt(0), + fElectronMinPt(0), fPHOSMinPhi(219.), fPHOSMaxPhi(321.), fPHOSEta(0.13), @@ -133,7 +142,6 @@ AliGenPythia::AliGenPythia(): { // Default Constructor fEnergyCMS = 5500.; - SetNuclei(0,0); if (!AliPythiaRndm::GetPythiaRandom()) AliPythiaRndm::SetPythiaRandom(GetRandom()); } @@ -141,6 +149,7 @@ AliGenPythia::AliGenPythia(): AliGenPythia::AliGenPythia(Int_t npart) :AliGenMC(npart), fProcess(kPyCharm), + fItune(-1), fStrucFunc(kCTEQ5L), fKineBias(0.), fTrials(0), @@ -164,9 +173,13 @@ AliGenPythia::AliGenPythia(Int_t npart) fGinit(kTRUE), fGfinal(kTRUE), fHadronisation(kTRUE), + fPatchOmegaDalitz(0), fNpartons(0), fReadFromFile(kFALSE), fQuench(kFALSE), + fQhat(0.), + fLength(0.), + fImpact(0.), fPtKick(1.), fFullEvent(kTRUE), fDecayer(new AliDecayerPythia()), @@ -201,18 +214,21 @@ AliGenPythia::AliGenPythia(Int_t npart) fTriggerEta(0.9), fTriggerMultiplicity(0), fTriggerMultiplicityEta(0), + fTriggerMultiplicityPtMin(0), fCountMode(kCountAll), fHeader(0), fRL(0), - fFileName(0), + fkFileName(0), fFragPhotonInCalo(kFALSE), fPi0InCalo(kFALSE) , fPhotonInCalo(kFALSE), + fEleInEMCAL(kFALSE), fCheckEMCAL(kFALSE), fCheckPHOS(kFALSE), fCheckPHOSeta(kFALSE), fFragPhotonOrPi0MinPt(0), fPhotonMinPt(0), + fElectronMinPt(0), fPHOSMinPhi(219.), fPHOSMaxPhi(321.), fPHOSEta(0.13), @@ -231,7 +247,6 @@ AliGenPythia::AliGenPythia(Int_t npart) // Set random number generator if (!AliPythiaRndm::GetPythiaRandom()) AliPythiaRndm::SetPythiaRandom(GetRandom()); - SetNuclei(0,0); } AliGenPythia::~AliGenPythia() @@ -355,14 +370,14 @@ void AliGenPythia::Init() if (fReadFromFile) { - fRL = AliRunLoader::Open(fFileName, "Partons"); + fRL = AliRunLoader::Open(fkFileName, "Partons"); fRL->LoadKinematics(); fRL->LoadHeader(); } else { 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 @@ -418,7 +433,12 @@ void AliGenPythia::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: @@ -447,6 +467,7 @@ void AliGenPythia::Init() fParentSelect[0] = 443; break; case kPyMbDefault: + case kPyMbAtlasTuneMC09: case kPyMb: case kPyMbWithDirectPhoton: case kPyMbNonDiffr: @@ -495,10 +516,14 @@ void AliGenPythia::Init() Warning("Init","SetNuclei used. Use SetProjectile + SetTarget instead. fDyBoost has been reset to 0\n"); } - if (fQuench) { + fPythia->SetPARJ(200, 0.0); + fPythia->SetPARJ(199, 0.0); + fPythia->SetPARJ(198, 0.0); + fPythia->SetPARJ(197, 0.0); + + if (fQuench == 1) { fPythia->InitQuenching(0., 0.1, 0.6e6, 0); } - fPythia->SetPARJ(200, 0.0); if (fQuench == 3) { // Nestor's change of the splittings @@ -509,6 +534,19 @@ void AliGenPythia::Init() fPythia->SetMSTJ(47, 0); // No correction back to hard scattering element fPythia->SetMSTJ(50, 0); // No coherence in first branching fPythia->SetPARJ(82, 1.); // Cut off for parton showers + } else if (fQuench == 4) { + // Armesto-Cunqueiro-Salgado change of the splittings. + AliFastGlauber* glauber = AliFastGlauber::Instance(); + glauber->Init(2); + //read and store transverse almonds corresponding to differnt + //impact parameters. + glauber->SetCentralityClass(0.,0.1); + fPythia->SetPARJ(200, 1.); + fPythia->SetPARJ(198, fQhat); + fPythia->SetPARJ(199, fLength); + fPythia->SetMSTJ(42, 2); // angular ordering + fPythia->SetMSTJ(44, 2); // option to run alpha_s + fPythia->SetPARJ(82, 1.); // Cut off for parton showers } } @@ -545,6 +583,15 @@ void AliGenPythia::Generate() // Switch hadronisation off // fPythia->SetMSTJ(1, 0); + + if (fQuench ==4){ + Double_t bimp; + // Quenching comes through medium-modified splitting functions. + AliFastGlauber::Instance()->GetRandomBHard(bimp); + fPythia->SetPARJ(197, bimp); + fImpact = bimp; + fPythia->Qpygin0(); + } // // Either produce new event or read partons from file // @@ -582,11 +629,21 @@ void AliGenPythia::Generate() // // .. and perform hadronisation // printf("Calling hadronisation %d\n", fPythia->GetN()); - fPythia->Pyexec(); + + if (fPatchOmegaDalitz) { + fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 0); + fPythia->Pyexec(); + fPythia->DalitzDecays(); + fPythia->SetMDCY(fPythia->Pycomp(111) ,1, 1); + } + fPythia->Pyexec(); } fTrials++; fPythia->ImportParticles(&fParticles,"All"); - Boost(); + + if (TMath::Abs(fDyBoost) > 1.e-4) Boost(); + if(TMath::Abs(fXingAngleY) > 1.e-10) BeamCrossAngle(); + // // // @@ -613,6 +670,7 @@ void AliGenPythia::Generate() Int_t nTkbles = 0; // Trackable particles if (fProcess != kPyMbDefault && fProcess != kPyMb && + fProcess != kPyMbAtlasTuneMC09 && fProcess != kPyMbWithDirectPhoton && fProcess != kPyJets && fProcess != kPyDirectGamma && @@ -621,7 +679,8 @@ void AliGenPythia::Generate() fProcess != kPyW && fProcess != kPyZ && fProcess != kPyCharmppMNRwmi && - fProcess != kPyBeautyppMNRwmi) { + fProcess != kPyBeautyppMNRwmi && + fProcess != kPyBeautyJets) { for (i = 0; i < np; i++) { TParticle* iparticle = (TParticle *) fParticles.At(i); @@ -850,7 +909,7 @@ Int_t AliGenPythia::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 == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) { TParticle* jet1 = (TParticle *) fParticles.At(6); TParticle* jet2 = (TParticle *) fParticles.At(7); if (!CheckTrigger(jet1, jet2)) { @@ -866,7 +925,7 @@ Int_t AliGenPythia::GenerateMB() Int_t pdg = 0; if (fFragPhotonInCalo) pdg = 22 ; // Photon - else if (fPi0InCalo) pdg = 111 ; // Pi0 + else if (fPi0InCalo) pdg = 111 ; // Pi0 for (i=0; i< np; i++) { TParticle* iparticle = (TParticle *) fParticles.At(i); @@ -885,10 +944,40 @@ Int_t AliGenPythia::GenerateMB() } } } - if(!ok) + if(!ok) { + delete[] pParent; + return 0; + } + } + + // Select beauty jets with electron in EMCAL + if (fProcess == kPyBeautyJets && fEleInEMCAL) { + + Bool_t ok = kFALSE; + + Int_t pdg = 11; //electron + + Float_t pt = 0.; + Float_t eta = 0.; + Float_t phi = 0.; + for (i=0; i< np; i++) { + TParticle* iparticle = (TParticle *) fParticles.At(i); + if(iparticle->GetStatusCode()==1 && TMath::Abs(iparticle->GetPdgCode())==pdg && + iparticle->Pt() > fElectronMinPt){ + pt = iparticle->Pt(); + phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees + eta =TMath::Abs(iparticle->Eta()); //in calos etamin=-etamax + if(IsInEMCAL(phi,eta)) + ok =kTRUE; + } + } + if(!ok) { + delete[] pParent; return 0; + } + + AliDebug(5,Form("Found an electron jet (pt,eta,phi) = (%f,%f,%f)",pt,eta,phi)); } - // Check for minimum multiplicity if (fTriggerMultiplicity > 0) { Int_t multiplicity = 0; @@ -898,17 +987,15 @@ Int_t AliGenPythia::GenerateMB() Int_t statusCode = iparticle->GetStatusCode(); // Initial state particle - if (statusCode > 20) + if (statusCode != 1) continue; - - // skip quarks and gluons - Int_t pdgCode = TMath::Abs(iparticle->GetPdgCode()); - if (pdgCode <= 10 || pdgCode == 21) - continue; - + // eta cut if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta) continue; - + // pt cut + if (iparticle->Pt() < fTriggerMultiplicityPtMin) + continue; + TParticlePDG* pdgPart = iparticle->GetPDG(); if (pdgPart && pdgPart->Charge() == 0) continue; @@ -920,8 +1007,7 @@ Int_t AliGenPythia::GenerateMB() delete [] pParent; return 0; } - - Printf("Triggered on event with multiplicity of %d > %d", multiplicity, fTriggerMultiplicity); + Printf("Triggered on event with multiplicity of %d >= %d", multiplicity, fTriggerMultiplicity); } // Select events with a photon pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi @@ -954,8 +1040,10 @@ Int_t AliGenPythia::GenerateMB() if(!okd && iphcand != -1) // execute rotation in phi RotatePhi(iphcand,okd); - if(!okd) - return 0; + if(!okd) { + delete [] pParent; + return 0; + } } if (fTriggerParticle) { @@ -978,7 +1066,8 @@ Int_t AliGenPythia::GenerateMB() // 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) { + + if (fProcess == kPyCharmppMNRwmi || fProcess == kPyBeautyppMNRwmi || fProcess == kPyBeautyJets) { TParticle *partCheck; TParticle *mother; Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE; @@ -1024,6 +1113,7 @@ Int_t AliGenPythia::GenerateMB() fProcess == kPyZ || fProcess == kPyMbDefault || fProcess == kPyMb || + fProcess == kPyMbAtlasTuneMC09 || fProcess == kPyMbWithDirectPhoton || fProcess == kPyMbNonDiffr) && (fCutOnChild == 1) ) { @@ -1042,7 +1132,7 @@ Int_t AliGenPythia::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 == kPyBeautyJets) && ks == 21 && km == 0 && i>1)) { nc++; if (ks == 1) trackIt = 1; Int_t ipa = iparticle->GetFirstMother()-1; @@ -1144,14 +1234,15 @@ void AliGenPythia::MakeHeader() // // Event Vertex fHeader->SetPrimaryVertex(fVertex); - + fHeader->SetInteractionTime(fEventTime); // // Number of primaries fHeader->SetNProduced(fNprimaries); // // Jets that have triggered - if (fProcess == kPyJets || fProcess == kPyDirectGamma) + //Need to store jets for b-jet studies too! + if (fProcess == kPyJets || fProcess == kPyDirectGamma || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) { Int_t ntrig, njet; Float_t jets[4][10]; @@ -1183,12 +1274,14 @@ void AliGenPythia::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); - } else { + } else if (fQuench == 2){ // Pyquen Double_t r1 = PARIMP.rb1; Double_t r2 = PARIMP.rb2; @@ -1198,10 +1291,20 @@ void AliGenPythia::MakeHeader() xp = r * TMath::Cos(phi); yp = r * TMath::Sin(phi); + } else if (fQuench == 4) { + // QPythia + Double_t xy[2]; + Double_t i0i1[2]; + AliFastGlauber::Instance()->GetSavedXY(xy); + AliFastGlauber::Instance()->GetSavedI0I1(i0i1); + xp = xy[0]; + yp = xy[1]; + ((AliGenPythiaEventHeader*) fHeader)->SetImpactParameter(fImpact); } + ((AliGenPythiaEventHeader*) fHeader)->SetXYJet(xp, yp); ((AliGenPythiaEventHeader*) fHeader)->SetZQuench(z); - } + } // // Store pt^hard ((AliGenPythiaEventHeader*) fHeader)->SetPtHard(fPythia->GetVINT(47)); @@ -1227,7 +1330,7 @@ Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2) pdg[1] = jet2->GetPdgCode(); Bool_t triggered = kFALSE; - if (fProcess == kPyJets) { + if (fProcess == kPyJets || fProcess == kPyBeautyJets || fProcess == kPyBeautyppMNRwmi) { Int_t njets = 0; Int_t ntrig = 0; Float_t jets[4][10]; @@ -1364,17 +1467,19 @@ void AliGenPythia::LoadEvent(TObjArray* stack, Int_t flag, Int_t reHadr) for (Int_t part = 0; part < npart; part++) { TParticle *mPart = dynamic_cast(stack->At(part)); + if (!mPart) continue; + 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) { - ks -= 10; - idf = -1; - idl = -1; - } + if (ks == 11 || ks == 12) { + ks -= 10; + idf = -1; + idl = -1; + } } Float_t px = mPart->Px(); @@ -1537,8 +1642,8 @@ void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd) //loop for all particles and produce the phi rotation Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons; Double_t oldphi, newphi; - Double_t newVx, newVy, R, Vz, time; - Double_t newPx, newPy, pt, Pz, e; + Double_t newVx, newVy, r, vZ, time; + Double_t newPx, newPy, pt, pz, e; for(Int_t i=0; i< np; i++) { TParticle* iparticle = (TParticle *) fParticles.At(i); oldphi = iparticle->Phi(); @@ -1546,21 +1651,21 @@ void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd) if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle - R = iparticle->R(); - newVx = R*TMath::Cos(newphi); - newVy = R*TMath::Sin(newphi); - Vz = iparticle->Vz(); // don't transform + r = iparticle->R(); + newVx = r * TMath::Cos(newphi); + newVy = r * TMath::Sin(newphi); + vZ = iparticle->Vz(); // don't transform time = iparticle->T(); // don't transform pt = iparticle->Pt(); - newPx = pt*TMath::Cos(newphi); - newPy = pt*TMath::Sin(newphi); - Pz = iparticle->Pz(); // don't transform - e = iparticle->Energy(); // don't transform + newPx = pt * TMath::Cos(newphi); + newPy = pt * TMath::Sin(newphi); + pz = iparticle->Pz(); // don't transform + e = iparticle->Energy(); // don't transform // apply rotation - iparticle->SetProductionVertex(newVx, newVy, Vz, time); - iparticle->SetMomentum(newPx, newPy, Pz, e); + iparticle->SetProductionVertex(newVx, newVy, vZ, time); + iparticle->SetMomentum(newPx, newPy, pz, e); } //end particle loop @@ -1572,46 +1677,4 @@ void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd) } -#ifdef never -void AliGenPythia::Streamer(TBuffer &R__b) -{ - // Stream an object of class AliGenPythia. - - if (R__b.IsReading()) { - Version_t R__v = R__b.ReadVersion(); if (R__v) { } - AliGenerator::Streamer(R__b); - R__b >> (Int_t&)fProcess; - R__b >> (Int_t&)fStrucFunc; - R__b >> (Int_t&)fForceDecay; - R__b >> fEnergyCMS; - R__b >> fKineBias; - R__b >> fTrials; - fParentSelect.Streamer(R__b); - fChildSelect.Streamer(R__b); - R__b >> fXsection; -// (AliPythia::Instance())->Streamer(R__b); - R__b >> fPtHardMin; - R__b >> fPtHardMax; -// if (fDecayer) fDecayer->Streamer(R__b); - } else { - R__b.WriteVersion(AliGenPythia::IsA()); - AliGenerator::Streamer(R__b); - R__b << (Int_t)fProcess; - R__b << (Int_t)fStrucFunc; - R__b << (Int_t)fForceDecay; - R__b << fEnergyCMS; - R__b << fKineBias; - R__b << fTrials; - fParentSelect.Streamer(R__b); - fChildSelect.Streamer(R__b); - R__b << fXsection; -// R__b << fPythia; - R__b << fPtHardMin; - R__b << fPtHardMax; - // fDecayer->Streamer(R__b); - } -} -#endif - -