#include "AliStack.h"
#include "AliRunLoader.h"
#include "AliMC.h"
+#include "AliLog.h"
#include "PyquenCommon.h"
+#include "AliLog.h"
ClassImp(AliGenPythia)
fPHOSEta(0.13),
fEMCALMinPhi(79.),
fEMCALMaxPhi(191.),
- fEMCALEta(0.71)
-
+ fEMCALEta(0.71),
+ fkTuneForDiff(0),
+ fProcDiff(0)
{
// Default Constructor
fEnergyCMS = 5500.;
fPHOSEta(0.13),
fEMCALMinPhi(79.),
fEMCALMaxPhi(191.),
- fEMCALEta(0.71)
+ fEMCALEta(0.71),
+ fkTuneForDiff(0),
+ fProcDiff(0)
{
// default charm production at 5. 5 TeV
// semimuonic decay
if (!fPythia) fPythia=(AliPythia*) fMCEvGen;
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;
//
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 = iparticle->GetFirstMother()-1;
Int_t iparent = (ipa > -1) ? pParent[ipa] : -1;
//
Int_t i, kf, nt, iparent;
Int_t nc = 0;
- Float_t p[4];
- Float_t polar[3] = {0,0,0};
- Float_t origin[3] = {0,0,0};
+ Double_t p[4];
+ Double_t polar[3] = {0,0,0};
+ Double_t origin[3] = {0,0,0};
// converts from mm/c to s
const Float_t kconv=0.001/2.999792458e8;
AliDebug(5,Form("Found an electron jet (pt,eta,phi) = (%f,%f,%f)",pt,eta,phi));
}
+ // Check for diffraction
+ if(fkTuneForDiff) {
+ if(fItune==320) {
+ if(!CheckDiffraction()) {
+ delete [] pParent;
+ return 0;
+ }
+ }
+ }
+
// Check for minimum multiplicity
if (fTriggerMultiplicity > 0) {
Int_t multiplicity = 0;
if(!okd && iphcand != -1) // execute rotation in phi
RotatePhi(iphcand,okd);
- if(!okd)
- return 0;
+ if(!okd) {
+ delete [] pParent;
+ return 0;
+ }
}
if (fTriggerParticle) {
TParticle *mother;
Bool_t theQ=kFALSE,theQbar=kFALSE,inYcut=kFALSE;
Bool_t theChild=kFALSE;
+ Bool_t triggered=kFALSE;
Float_t y;
Int_t pdg,mpdg,mpdgUpperFamily;
for(i=0; i<np; i++) {
(partCheck->Energy()-partCheck->Pz()+1.e-13));
if(y>fYMin && y<fYMax) inYcut=kTRUE;
}
+ if(fTriggerParticle) {
+ if(TMath::Abs(pdg)==fTriggerParticle) triggered=kTRUE;
+ }
if(fCutOnChild && TMath::Abs(pdg) == fPdgCodeParticleforAcceptanceCut) {
Int_t mi = partCheck->GetFirstMother() - 1;
if(mi<0) continue;
delete[] pParent;
return 0;
}
+ if(fTriggerParticle && !triggered) { // particle requested is not produced
+ delete[] pParent;
+ return 0;
+ }
}
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],
fHeader = new AliGenPythiaEventHeader("Pythia");
//
// Event type
+ if(fProcDiff>0){
+ // if(fProcDiff == 92 || fProcDiff == 93) printf("\n\n\n\n\n");
+ // printf("fPythia->GetMSTI(1) = %d fProcDiff = %d\n",fPythia->GetMSTI(1), fProcDiff);
+ ((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fProcDiff);
+ }
+ else
((AliGenPythiaEventHeader*) fHeader)->SetProcessType(fPythia->GetMSTI(1));
//
// Number of trials
//
// Event Vertex
fHeader->SetPrimaryVertex(fVertex);
- fHeader->SetInteractionTime(fEventTime);
+ fHeader->SetInteractionTime(fTime+fEventTime);
//
// Number of primaries
fHeader->SetNProduced(fNprimaries);
fHeader = 0x0;
}
-Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
+Bool_t AliGenPythia::CheckTrigger(const TParticle* jet1, const TParticle* jet2)
{
// Check the kinematic trigger condition
//
}
}
-void AliGenPythia::LoadEvent(TObjArray* stack, Int_t flag, Int_t reHadr)
+void AliGenPythia::LoadEvent(const TObjArray* stack, Int_t flag, Int_t reHadr)
{
//
// Load event into Pythia Common Block
return;
}
-Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta)
+Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta) const
{
// Is particle in EMCAL acceptance?
// phi in degrees, etamin=-etamax
return kFALSE;
}
-Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta)
+Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta) const
{
// Is particle in PHOS acceptance?
// Acceptance slightly larger considered.
//calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
- Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
+ Double_t phiPHOS = (AliPythiaRndm::GetPythiaRandom())->Uniform(phiPHOSmin,phiPHOSmax);
//calculate deltaphi
TParticle* ph = (TParticle *) fParticles.At(iphcand);
+Bool_t AliGenPythia::CheckDiffraction()
+{
+ // use this method only with Perugia-0 tune!
+
+ // printf("AAA\n");
+
+ Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
+
+ Int_t iPart1=-1;
+ Int_t iPart2=-1;
+
+ Double_t y1 = 1e10;
+ Double_t y2 = -1e10;
+
+ const Int_t kNstable=20;
+ const Int_t pdgStable[20] = {
+ 22, // Photon
+ 11, // Electron
+ 12, // Electron Neutrino
+ 13, // Muon
+ 14, // Muon Neutrino
+ 15, // Tau
+ 16, // Tau Neutrino
+ 211, // Pion
+ 321, // Kaon
+ 311, // K0
+ 130, // K0s
+ 310, // K0l
+ 2212, // Proton
+ 2112, // Neutron
+ 3122, // Lambda_0
+ 3112, // Sigma Minus
+ 3222, // Sigma Plus
+ 3312, // Xsi Minus
+ 3322, // Xsi0
+ 3334 // Omega
+ };
+
+ for (Int_t i = 0; i < np; i++) {
+ TParticle * part = (TParticle *) fParticles.At(i);
+
+ Int_t statusCode = part->GetStatusCode();
+
+ // Initial state particle
+ if (statusCode != 1)
+ continue;
+
+ Int_t pdg = TMath::Abs(part->GetPdgCode());
+ Bool_t isStable = kFALSE;
+ for (Int_t i1 = 0; i1 < kNstable; i1++) {
+ if (pdg == pdgStable[i1]) {
+ isStable = kTRUE;
+ break;
+ }
+ }
+ if(!isStable)
+ continue;
+
+ Double_t y = part->Y();
+
+ if (y < y1)
+ {
+ y1 = y;
+ iPart1 = i;
+ }
+ if (y > y2)
+ {
+ y2 = y;
+ iPart2 = i;
+ }
+ }
+
+ if(iPart1<0 || iPart2<0) return kFALSE;
+
+ y1=TMath::Abs(y1);
+ y2=TMath::Abs(y2);
+
+ TParticle * part1 = (TParticle *) fParticles.At(iPart1);
+ TParticle * part2 = (TParticle *) fParticles.At(iPart2);
+
+ Int_t pdg1 = part1->GetPdgCode();
+ Int_t pdg2 = part2->GetPdgCode();
+
+
+ Int_t iPart = -1;
+ if (pdg1 == 2212 && pdg2 == 2212)
+ {
+ if(y1 > y2)
+ iPart = iPart1;
+ else if(y1 < y2)
+ iPart = iPart2;
+ else {
+ iPart = iPart1;
+ if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)>0.5) iPart = iPart2;
+ }
+ }
+ else if (pdg1 == 2212)
+ iPart = iPart1;
+ else if (pdg2 == 2212)
+ iPart = iPart2;
+
+
+
+
+
+ Double_t M=-1.;
+ if(iPart>0) {
+ TParticle * part = (TParticle *) fParticles.At(iPart);
+ Double_t E= part->Energy();
+ Double_t P= part->P();
+ M= TMath::Sqrt((fEnergyCMS-E-P)*(fEnergyCMS-E+P));
+ }
+
+const Int_t nbin=120;
+Double_t bin[]={
+1.080000, 1.274258, 1.468516, 1.662773, 1.857031, 2.051289,
+2.245547, 2.439805, 2.634062, 2.828320, 3.022578, 3.216836,
+3.411094, 3.605352, 3.799609, 3.993867, 4.188125, 4.382383,
+4.576641, 4.770898, 4.965156, 5.547930, 6.130703, 6.713477,
+7.296250, 7.879023, 8.461797, 9.044570, 9.627344, 10.210117,
+10.792891, 11.375664, 11.958437, 12.541211, 13.123984, 13.706758,
+14.289531, 14.872305, 15.455078, 16.037852, 16.620625, 17.203398,
+17.786172, 18.368945, 18.951719, 19.534492, 20.117266, 20.700039,
+21.282812, 21.865586, 22.448359, 23.031133, 23.613906, 24.196680,
+24.779453, 25.362227, 25.945000, 26.527773, 27.110547, 27.693320,
+28.276094, 28.858867, 29.441641, 30.024414, 30.607187, 31.189961,
+31.772734, 32.355508, 32.938281, 33.521055, 34.103828, 34.686602,
+35.269375, 35.852148, 36.434922, 37.017695, 37.600469, 38.183242,
+38.766016, 39.348789, 39.931562, 40.514336, 41.097109, 41.679883,
+42.262656, 42.845430, 43.428203, 44.010977, 44.593750, 45.176523,
+45.759297, 46.342070, 46.924844, 47.507617, 48.090391, 48.673164,
+49.255937, 49.838711, 50.421484, 57.220508, 64.019531, 70.818555,
+77.617578, 84.416602, 91.215625, 98.014648, 104.813672, 111.612695,
+118.411719, 125.210742, 132.009766, 138.808789, 145.607812, 152.406836,
+159.205859, 166.004883, 172.803906, 179.602930, 186.401953, 193.200977,
+200.000000};
+Double_t w[]={
+1.000000, 0.275812, 0.267308, 0.263268, 0.263828, 0.263039,
+0.260546, 0.259344, 0.255477, 0.253782, 0.253562, 0.252492,
+0.251076, 0.247862, 0.248933, 0.243599, 0.244255, 0.238506,
+0.239546, 0.237845, 0.237977, 0.229379, 0.224771, 0.217581,
+0.208860, 0.207241, 0.198955, 0.196449, 0.193827, 0.190220,
+0.184715, 0.183067, 0.178325, 0.175887, 0.171299, 0.168718,
+0.167453, 0.165153, 0.163457, 0.159637, 0.156855, 0.155488,
+0.154545, 0.155968, 0.150488, 0.148797, 0.145358, 0.146196,
+0.145891, 0.142752, 0.145072, 0.141265, 0.141857, 0.138462,
+0.142992, 0.141357, 0.139391, 0.139629, 0.135197, 0.135731,
+0.133194, 0.137190, 0.135745, 0.134522, 0.136094, 0.130405,
+0.128371, 0.131680, 0.130591, 0.133539, 0.129370, 0.128254,
+0.128262, 0.131088, 0.128294, 0.130070, 0.124553, 0.131489,
+0.128038, 0.122997, 0.130699, 0.125630, 0.124746, 0.123679,
+0.127864, 0.125776, 0.126272, 0.121492, 0.124929, 0.122221,
+0.127017, 0.123706, 0.122701, 0.123731, 0.122219, 0.120783,
+0.120324, 0.120228, 0.123615, 0.120589, 0.119549, 0.118836,
+0.118146, 0.120175, 0.122031, 0.122076, 0.122849, 0.122627,
+0.126281, 0.127696, 0.129084, 0.130001, 0.134062, 0.134786,
+0.137286, 0.136625, 0.139091, 0.143692, 0.144781, 0.146768};
+
+ Double_t wSD=1.;
+ Double_t wDD=0.178783;
+ //Double_t wND=0.220200;
+ Double_t wND=0.220200+0.001;
+
+ if(M>-1 && M<bin[0]) return kFALSE;
+ if(M>bin[nbin]) M=-1;
+
+ Int_t procType=fPythia->GetMSTI(1);
+ Int_t proc0=2;
+ if(procType== 94) proc0=1;
+ if(procType== 92 || procType== 93) proc0=0;
+
+
+ // printf("M = %f bin[nbin] = %f\n",M, bin[nbin]);
+
+ Int_t proc=2;
+ if(M>0) proc=0;
+ else if(proc0==1) proc=1;
+
+ if(proc==0 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wSD) return kFALSE;
+ if(proc==1 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wDD) return kFALSE;
+ if(proc==2 && (AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.) > wND) return kFALSE;
+
+
+ // if(proc==1 || proc==2) return kFALSE;
+
+ if(proc!=0) {
+ if(proc0!=0) fProcDiff = procType;
+ else fProcDiff = 95;
+ return kTRUE;
+ }
+
+ Int_t ibin=nbin-1;
+ for(Int_t i=1; i<=nbin; i++)
+ if(M<=bin[i]) {
+ ibin=i-1;
+ // printf("Mi> %f && Mi< %f\n", bin[i-1], bin[i]);
+ break;
+ }
+
+ // printf("w[ibin] = %f\n", w[ibin]);
+
+ if((AliPythiaRndm::GetPythiaRandom())->Uniform(0.,1.)> w[ibin]) return kFALSE;
+
+ // printf("iPart = %d\n", iPart);
+
+ if(iPart==iPart1) fProcDiff=93;
+ else if(iPart==iPart2) fProcDiff=92;
+ else {
+ printf("EROOR: iPart!=iPart1 && iPart!=iPart2\n");
+
+ }
+
+ return kTRUE;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff