X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=TTherminator%2FAliGenTherminator.cxx;h=4f42251c48654f08098177703d9ac71f3702b2bc;hp=dffe99c7b6ff33ee26d1b0bad80ee26638361e91;hb=f2cef1b5b1d8646af95064ccd41b3bab8a368707;hpb=2e9679195eb6ff424774adcc2010861a8517fdaf diff --git a/TTherminator/AliGenTherminator.cxx b/TTherminator/AliGenTherminator.cxx index dffe99c7b6f..4f42251c486 100644 --- a/TTherminator/AliGenTherminator.cxx +++ b/TTherminator/AliGenTherminator.cxx @@ -1,6 +1,8 @@ // ALICE event generator based on the THERMINATOR model // It reads the test output of the model and puts it onto // the stack +// It has an option to use the Lhyquid3D input freeze-out +// hypersurface // Author: Adam.Kisiel@cern.ch #include @@ -14,6 +16,7 @@ #include "AliConst.h" #include "AliDecayer.h" #include "AliGenEventHeader.h" +#include "AliGenHijingEventHeader.h" #include "AliGenTherminator.h" #include "AliLog.h" #include "AliRun.h" @@ -28,6 +31,7 @@ AliGenTherminator::AliGenTherminator(): fEventNumber(0), fFileName(""), fFreezeOutModel(""), + fFOHSlocation(""), fTemperature(0.1656), fMiuI(-0.0009), fMiuS(0.0), @@ -41,7 +45,15 @@ AliGenTherminator::AliGenTherminator(): fBWDelay(0.0) { // Default constructor - fParticles = new TClonesArray("TParticle",20000); + fFOHSlocation = ""; + + fEnergyCMS = 5500.; + fAProjectile = 208; + fZProjectile = 82; + fProjectile = "A"; + fATarget = 208; + fZTarget = 82; + fTarget = "A"; } AliGenTherminator::AliGenTherminator(Int_t npart): AliGenMC(npart), @@ -49,6 +61,7 @@ AliGenTherminator::AliGenTherminator(Int_t npart): fEventNumber(0), fFileName(""), fFreezeOutModel(""), + fFOHSlocation(""), fTemperature(0.1656), fMiuI(-0.0009), fMiuS(0.0), @@ -62,8 +75,14 @@ AliGenTherminator::AliGenTherminator(Int_t npart): fBWDelay(0.0) { // Constructor specifying the size of the particle table - fParticles = new TClonesArray("TParticle",20000); fNprimaries = 0; + fEnergyCMS = 5500.; + fAProjectile = 208; + fZProjectile = 82; + fProjectile = "A"; + fATarget = 208; + fZTarget = 82; + fTarget = "A"; } AliGenTherminator::~AliGenTherminator() @@ -94,48 +113,131 @@ void AliGenTherminator::Generate() ((TTherminator *) fMCEvGen)->GenerateEvent(); AliWarning("Generated"); - ((TTherminator *) fMCEvGen)->ImportParticles(fParticles); + ((TTherminator *) fMCEvGen)->ImportParticles(&fParticles); - Int_t np = fParticles->GetEntriesFast(); + Int_t np = fParticles.GetEntriesFast(); AliWarning(Form("Imported %d particles", np)); + Int_t *idsOnStack; + idsOnStack = new Int_t[np]; + TParticle *iparticle; + Double_t evrot = gRandom->Rndm()*TMath::Pi(); for (int i = 0; i < np; i++) { - iparticle = (TParticle *) fParticles->At(i); + iparticle = (TParticle *) fParticles.At(i); Bool_t hasMother = (iparticle->GetFirstMother() >=0); Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0); - kf = iparticle->GetPdgCode(); - ks = iparticle->GetStatusCode(); - p[0] = iparticle->Px(); - p[1] = iparticle->Py(); - p[2] = iparticle->Pz(); - mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass(); - energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]); - origin[0] = origin0[0]+iparticle->Vx(); - origin[1] = origin0[1]+iparticle->Vy(); - origin[2] = origin0[2]+iparticle->Vz(); - - imo = -1; - TParticle* mother = 0; - if (hasMother) { - imo = iparticle->GetFirstMother(); - mother = (TParticle *) fParticles->At(imo); - } // if has mother - Bool_t tFlag = (!hasDaughter); - - printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo); - PushTrack(tFlag,imo,kf, - p[0],p[1],p[2],energy, - origin[0],origin[1],origin[2],iparticle->T()*0.197327*1e-13/300000000, - polar[0],polar[1],polar[2], - hasMother ? kPDecay:kPNoProcess,nt); - - fNprimaries++; - KeepTrack(nt); + if (hasDaughter) { + // This particle has decayed + // It will not be tracked + // Add it only once with coorduinates not + // smeared with primary vertex position + + kf = iparticle->GetPdgCode(); + ks = iparticle->GetStatusCode(); + Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px()); + Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py()); + p[0] = arho*TMath::Cos(aphi + evrot); + p[1] = arho*TMath::Sin(aphi + evrot); + // p[0] = iparticle->Px(); + // p[1] = iparticle->Py(); + p[2] = iparticle->Pz(); + mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass(); + energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]); + + Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx()); + Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy()); + origin[0] = vrho*TMath::Cos(vphi + evrot); + origin[1] = vrho*TMath::Sin(vphi + evrot); + origin[2] = iparticle->Vz(); + + imo = -1; + // TParticle* mother = 0; + if (hasMother) { + imo = iparticle->GetFirstMother(); + // mother = (TParticle *) fParticles.At(imo); + } // if has mother + Bool_t tFlag = (!hasDaughter); + + // printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo>=0?idsOnStack[imo]:imo); + PushTrack(tFlag,imo>=0?idsOnStack[imo]:imo,kf, + p[0],p[1],p[2],energy, + origin[0],origin[1],origin[2],iparticle->T(), + polar[0],polar[1],polar[2], + hasMother ? kPDecay:kPNoProcess,nt); + idsOnStack[i] = nt; + fNprimaries++; + KeepTrack(nt); + } + else { + // This is a final state particle + // It will be tracked + // Add it TWICE to the stack !!! + // First time with event-wide coordicates (for femtoscopy) - + // this one will not be tracked + // Second time with event-wide ccordiantes and vertex smearing + // this one will be tracked + + kf = iparticle->GetPdgCode(); + ks = iparticle->GetStatusCode(); + Double_t aphi = TMath::ATan2(iparticle->Py(), iparticle->Px()); + Double_t arho = TMath::Hypot(iparticle->Px(), iparticle->Py()); + p[0] = arho*TMath::Cos(aphi + evrot); + p[1] = arho*TMath::Sin(aphi + evrot); + // p[0] = iparticle->Px(); + // p[1] = iparticle->Py(); + p[2] = iparticle->Pz(); + mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass(); + energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]); + + Double_t vphi = TMath::ATan2(iparticle->Vy(), iparticle->Vx()); + Double_t vrho = TMath::Hypot(iparticle->Vx(), iparticle->Vy()); + origin[0] = vrho*TMath::Cos(vphi + evrot); + origin[1] = vrho*TMath::Sin(vphi + evrot); + origin[2] = iparticle->Vz(); + + imo = -1; + // TParticle* mother = 0; + if (hasMother) { + imo = iparticle->GetFirstMother(); + // mother = (TParticle *) fParticles.At(imo); + } // if has mother + Bool_t tFlag = (hasDaughter); + +// printf("Found mother %i with true id %i\n", imo, imo>=0?idsOnStack[imo]:imo); +// printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo>=0?idsOnStack[imo]:imo); + PushTrack(tFlag,imo>=0?idsOnStack[imo]:imo,kf, + p[0],p[1],p[2],energy, + origin[0],origin[1],origin[2],iparticle->T(), + polar[0],polar[1],polar[2], + hasMother ? kPDecay:kPNoProcess,nt); + idsOnStack[i] = nt; + fNprimaries++; + KeepTrack(nt); + + origin[0] = origin0[0]+vrho*TMath::Cos(vphi + evrot); + origin[1] = origin0[1]+vrho*TMath::Sin(vphi + evrot); + origin[2] = origin0[2]+iparticle->Vz(); + + imo = nt; + // mother = (TParticle *) fParticles.At(nt); + tFlag = (!hasDaughter); + +// printf("Pushing Track %d with status %d mother %d\n", kf, tFlag, imo); + PushTrack(tFlag,imo,kf, + p[0],p[1],p[2],energy, + origin[0],origin[1],origin[2],iparticle->T(), + polar[0],polar[1],polar[2], + hasMother ? kPDecay:kPNoProcess,nt); + fNprimaries++; + KeepTrack(nt); + } } + + SetHighWaterMark(fNprimaries); TArrayF eventVertex; @@ -144,12 +246,59 @@ void AliGenTherminator::Generate() eventVertex[1] = origin0[1]; eventVertex[2] = origin0[2]; +// Builds the event header, to be called after each event + AliGenEventHeader* header = new AliGenHijingEventHeader("Therminator"); + // Header - AliGenEventHeader* header = new AliGenEventHeader("Therminator"); + // AliGenEventHeader* header = new AliGenEventHeader("Therminator"); // Event Vertex - header->SetPrimaryVertex(eventVertex); - header->SetNProduced(fNprimaries); - gAlice->SetGenEventHeader(header); +// header->SetPrimaryVertex(eventVertex); +// header->SetNProduced(fNprimaries); + + ((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries); + ((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex); + ((AliGenHijingEventHeader*) header)->SetImpactParameter(0.0); + ((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0); + ((AliGenHijingEventHeader*) header)->SetHardScatters(0); + ((AliGenHijingEventHeader*) header)->SetParticipants(0, 0); + ((AliGenHijingEventHeader*) header)->SetCollisions(0, 0, 0, 0); + ((AliGenHijingEventHeader*) header)->SetSpectators(0, 0, 0, 0); + ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(evrot); + + +// 4-momentum vectors of the triggered jets. +// +// Before final state gluon radiation. +// TLorentzVector* jet1 = new TLorentzVector(fHijing->GetHINT1(21), +// fHijing->GetHINT1(22), +// fHijing->GetHINT1(23), +// fHijing->GetHINT1(24)); + +// TLorentzVector* jet2 = new TLorentzVector(fHijing->GetHINT1(31), +// fHijing->GetHINT1(32), +// fHijing->GetHINT1(33), +// fHijing->GetHINT1(34)); +// // After final state gluon radiation. +// TLorentzVector* jet3 = new TLorentzVector(fHijing->GetHINT1(26), +// fHijing->GetHINT1(27), +// fHijing->GetHINT1(28), +// fHijing->GetHINT1(29)); + +// TLorentzVector* jet4 = new TLorentzVector(fHijing->GetHINT1(36), +// fHijing->GetHINT1(37), +// fHijing->GetHINT1(38), +// fHijing->GetHINT1(39)); +// ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4); +// Bookkeeping for kinematic bias +// ((AliGenHijingEventHeader*) header)->SetTrials(fTrials); +// Event Vertex + header->SetPrimaryVertex(fVertex); + AddHeader(header); + fCollisionGeometry = (AliGenHijingEventHeader*) header; + + delete [] idsOnStack; + + // gAlice->SetGenEventHeader(header); } void AliGenTherminator::Init() @@ -188,7 +337,10 @@ void AliGenTherminator::ReadShareParticleTable() TParticlePDG *tParticleType; AliWarning(Form("Reading particle types from particles.data")); - ifstream in("particles.data"); + + TString aroot = gSystem->Getenv("ALICE_ROOT"); + ifstream in((aroot+"/TTherminator/data/SHARE/particles.data").Data()); + // ifstream in("particles.data"); int charge; @@ -235,12 +387,14 @@ void AliGenTherminator::ReadShareParticleTable() void AliGenTherminator::CreateTherminatorInputFile() { // Create Therminator input file + const char *aroot = gSystem->Getenv("ALICE_ROOT"); ofstream *ostr = new ofstream("therminator.in"); (*ostr) << "NumberOfEvents = 1" << endl; (*ostr) << "Randomize = 1" << endl; (*ostr) << "TableType = SHARE" << endl; - (*ostr) << "InputDirSHARE = ." << endl; + (*ostr) << "InputDirSHARE = "<< aroot << "/TTherminator/data/SHARE" << endl; (*ostr) << "EventOutputFile = " << fFileName.Data() << endl; + (*ostr) << "FOHSLocation = " << fFOHSlocation.Data() << endl; (*ostr) << "FreezeOutModel = " << fFreezeOutModel.Data() << endl; (*ostr) << "BWVt = " << fBWVt << endl; (*ostr) << "Tau = " << fTau << endl; @@ -258,4 +412,126 @@ void AliGenTherminator::SetModel(const char *model) { // Set the freeze-out model to use fFreezeOutModel = model; + AliWarning(Form("Selected model %s", fFreezeOutModel.Data())); + AliWarning(Form("FOHSLocation is %s", fFOHSlocation.Data())); +} + +void AliGenTherminator::SetLhyquidSet(const char *set) +{ + // Select one of pregenerated Lhyquid hypersurfaces + const char *aroot = gSystem->Getenv("ALICE_ROOT"); + if (strstr(set, "LHC500C0005")) { + AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface")); + AliWarning(Form(" Pb-Pb collisions, centrality 0-5 percent")); + AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV")); + AliWarning(Form(" freeze-out criteria Tf=145 MeV")); + AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C0005/FO.txt")); + fFOHSlocation.Append(aroot); + fFOHSlocation.Append("/TTherminator/data/LHC500C0005"); + } + if (strstr(set, "LHC500C0510")) { + AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface")); + AliWarning(Form(" Pb-Pb collisions, centrality 5-10 percent")); + AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV")); + AliWarning(Form(" freeze-out criteria Tf=145 MeV")); + AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C0510/FO.txt")); + fFOHSlocation.Append(aroot); + fFOHSlocation.Append("/TTherminator/data/LHC500C0510"); + } + if (strstr(set, "LHC500C1020")) { + AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface")); + AliWarning(Form(" Pb-Pb collisions, centrality 10-20 percent")); + AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV")); + AliWarning(Form(" freeze-out criteria Tf=145 MeV")); + AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C1020/FO.txt")); + fFOHSlocation.Append(aroot); + fFOHSlocation.Append("/TTherminator/data/LHC500C1020"); + } + else if (strstr(set, "LHC500C2030")) { + AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface")); + AliWarning(Form(" Pb-Pb collisions, centrality 20-30 percent")); + AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV")); + AliWarning(Form(" freeze-out criteria Tf=145 MeV")); + AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C2030/FO.txt")); + fFOHSlocation.Append(aroot); + fFOHSlocation.Append("/TTherminator/data/LHC500C2030"); + } + else if (strstr(set, "LHC500C3040")) { + AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface")); + AliWarning(Form(" Pb-Pb collisions, centrality 30-40 percent")); + AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV")); + AliWarning(Form(" freeze-out criteria Tf=145 MeV")); + AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C3040/FO.txt")); + fFOHSlocation.Append(aroot); + fFOHSlocation.Append("/TTherminator/data/LHC500C3040"); + } + else if (strstr(set, "LHC500C4050")) { + AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface")); + AliWarning(Form(" Pb-Pb collisions, centrality 40-50 percent")); + AliWarning(Form(" initial temperature at tau=1 fm in the center Ti=500 MeV")); + AliWarning(Form(" freeze-out criteria Tf=145 MeV")); + AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC500C4050/FO.txt")); + fFOHSlocation.Append(aroot); + fFOHSlocation.Append("/TTherminator/data/LHC500C4050"); + } + else if (strstr(set, "LHC276TC0005")) { + AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface")); + AliWarning(Form(" Pb-Pb collisions, 2.76TeV, centrality 00-05 percent")); + AliWarning(Form(" freeze-out criteria Tf=145 MeV")); + AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC0005/FO.txt")); + fFOHSlocation.Append(aroot); + fFOHSlocation.Append("/TTherminator/data/LHC276TC0005"); + } + else if (strstr(set, "LHC276TC0510")) { + AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface")); + AliWarning(Form(" Pb-Pb collisions, 2.76TeV, centrality 05-10 percent")); + AliWarning(Form(" freeze-out criteria Tf=145 MeV")); + AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC0510/FO.txt")); + fFOHSlocation.Append(aroot); + fFOHSlocation.Append("/TTherminator/data/LHC276TC0510"); + } + else if (strstr(set, "LHC276TC1020")) { + AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface")); + AliWarning(Form(" Pb-Pb collisions, 2.76TeV, centrality 10-20 percent")); + AliWarning(Form(" freeze-out criteria Tf=145 MeV")); + AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC1020/FO.txt")); + fFOHSlocation.Append(aroot); + fFOHSlocation.Append("/TTherminator/data/LHC276TC1020"); + } + else if (strstr(set, "LHC276TC2030")) { + AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface")); + AliWarning(Form(" Pb-Pb collisions, 2.76TeV, centrality 20-30 percent")); + AliWarning(Form(" freeze-out criteria Tf=145 MeV")); + AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC2030/FO.txt")); + fFOHSlocation.Append(aroot); + fFOHSlocation.Append("/TTherminator/data/LHC276TC2030"); + } + else if (strstr(set, "LHC276TC3040")) { + AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface")); + AliWarning(Form(" Pb-Pb collisions, 2.76TeV, centrality 30-40 percent")); + AliWarning(Form(" freeze-out criteria Tf=145 MeV")); + AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC3040/FO.txt")); + fFOHSlocation.Append(aroot); + fFOHSlocation.Append("/TTherminator/data/LHC276TC3040"); + } + else if (strstr(set, "LHC276TC4050")) { + AliWarning(Form("AliGenTherminator: Selected default Lhyquid hypersurface")); + AliWarning(Form(" Pb-Pb collisions, 2.76TeV, centrality 40-50 percent")); + AliWarning(Form(" freeze-out criteria Tf=145 MeV")); + AliWarning(Form(" for details see $(ALICE_ROOT)/TTherminator/data/LHC276TC4050/FO.txt")); + fFOHSlocation.Append(aroot); + fFOHSlocation.Append("/TTherminator/data/LHC276TC4050"); + } + else { + AliWarning(Form("Did not find Lhyquid set %s", set)); + AliWarning(Form("Reverting to default: current directory")); + fFOHSlocation += ""; + } + +} + +void AliGenTherminator::SetLhyquidInputDir(const char *inputdir) +{ + // Select Your own Lhyquid hypersurface + fFOHSlocation = inputdir; }