X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TUHKMgen%2FAliGenUHKM.cxx;h=8e19451dbab2934f8cf55793fc1999021ad77ba2;hb=46e461b837d309ea86cf11457ef75593250afcb3;hp=42a89ba2d6851630247e61a725497fc8140386c2;hpb=7b7936e9d4b5cb17c62a050675371a0242f22270;p=u%2Fmrichter%2FAliRoot.git diff --git a/TUHKMgen/AliGenUHKM.cxx b/TUHKMgen/AliGenUHKM.cxx index 42a89ba2d68..8e19451dbab 100755 --- a/TUHKMgen/AliGenUHKM.cxx +++ b/TUHKMgen/AliGenUHKM.cxx @@ -14,6 +14,9 @@ // // //////////////////////////////////////////////////////////////////////////// +#include +#include + #include "TUHKMgen.h" #ifndef DATABASE_PDG #include "DatabasePDG.h" @@ -26,8 +29,8 @@ #include #include #include -#include "TDatabasePDG.h" -#include "TSystem.h" +#include +#include #include "AliGenUHKM.h" #include "AliRun.h" @@ -37,15 +40,8 @@ #include "AliGenHijingEventHeader.h" #include "AliLog.h" -#include -#include using namespace std; -//class TClonesArray; -//enum TMCProcess; -//class AliRun; -//class TSystem; -//class AliDecayer; ClassImp(AliGenUHKM) @@ -59,7 +55,6 @@ AliGenUHKM::AliGenUHKM() { // Default constructor setting up default reasonable parameter values // for central Pb+Pb collisions at 5.5TeV - // cout << "AliGenUHKM::AliGenUHKM() IN" << endl; // LHC fHydjetParams.fSqrtS=5500; //LHC @@ -124,15 +119,14 @@ AliGenUHKM::AliGenUHKM() fHydjetParams.fIenglu=0; fHydjetParams.fIanglu=0; */ - strcpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT"))); - strcpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT"))); + strncpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")), 255); + strncpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")), 255); for(Int_t i=0; i<500; i++) { fStableFlagPDG[i] = 0; fStableFlagStatus[i] = kFALSE; } fStableFlagged = 0; - // cout << "AliGenUHKM::AliGenUHKM() OUT" << endl; } //_______________________________________ @@ -146,7 +140,7 @@ AliGenUHKM::AliGenUHKM(Int_t npart) // Constructor specifying the size of the particle table // and setting up default reasonable parameter values // for central Pb+Pb collisions at 5.5TeV - // cout << "AliGenUHKM::AliGenUHKM(Int_t) IN" << endl; + fName = "UHKM"; fTitle= "Particle Generator using UHKM 3.0"; @@ -215,21 +209,21 @@ AliGenUHKM::AliGenUHKM(Int_t npart) fHydjetParams.fIanglu=0; */ - strcpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT"))); - strcpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT"))); + strncpy(fParticleFilename, Form("%s/TUHKMgen/UHKM/particles.data", gSystem->Getenv("ALICE_ROOT")), 255); + strncpy(fDecayFilename, Form("%s/TUHKMgen/UHKM/tabledecay.txt", gSystem->Getenv("ALICE_ROOT")), 255); for(Int_t i=0; i<500; i++) { fStableFlagPDG[i] = 0; fStableFlagStatus[i] = kFALSE; } fStableFlagged = 0; - // cout << "AliGenUHKM::AliGenUHKM(Int_t) OUT" << endl; } //__________________________________________ AliGenUHKM::~AliGenUHKM() { // Destructor, do nothing + // delete fParticles; } void AliGenUHKM::SetAllParametersRHIC() @@ -310,19 +304,15 @@ void AliGenUHKM::Init() // further to the model. // HYDJET++ is initialized (average multiplicities are calculated, particle species definitions and decay // channels are loaded, etc.) - // cout << "AliGenUHKM::Init() IN" << endl; SetMC(new TUHKMgen()); fUHKMgen = (TUHKMgen*) fMCEvGen; SetAllParameters(); AliGenMC::Init(); - fUHKMgen->Initialize(); - CheckPDGTable(); // - - fUHKMgen->Print(); // - // cout << "AliGenUHKM::Init() OUT" << endl; + CheckPDGTable(); + fUHKMgen->Print(); } //________________________________________ @@ -330,98 +320,90 @@ void AliGenUHKM::Generate() { // Generate one HYDJET++ event, get the output and push particles further // to AliRoot's stack - // cout << "AliGenUHKM::Generate() IN" << endl; + Float_t polar[3] = {0,0,0}; Float_t origin[3] = {0,0,0}; Float_t origin0[3] = {0,0,0}; + Float_t time0 = 0.; + Float_t p[3]; Float_t v[3]; - Float_t mass, energy; + Float_t mass=0.0, energy=0.0; Vertex(); for(Int_t j=0; j<3; j++) origin0[j] = fVertex[j]; + time0 = fTime; - fTrials = 0; - // cout << "AliGenUHKM::Generate() fTrials = " << fTrials << endl; - - Int_t nt = 0; - + // Generate the event and import particles fUHKMgen->GenerateEvent(); - fTrials++; - fUHKMgen->ImportParticles(&fParticles,"All"); - Int_t np = fParticles.GetEntriesFast(); - // cout << "AliGenUHKM::Generate() GetEntries " <GetPdgCode(); - // if(kf==130 || kf==-130) - // cout << "AliGenUHKM::Generate() PDG = " << kf << endl; + // Generate a random phi used to rotate the whole event + Double_t eventRotation = gRandom->Rndm()*TMath::Pi(); + TParticle *iparticle; + Double_t partMomPhi=0.0, partPt=0.0; + Double_t partVtxPhi=0.0, partVtxR=0.0; + //_________ Loop for particles in the stack + for(Int_t i=0; iGetPdgCode(); Bool_t hasMother = (iparticle->GetFirstMother() >= 0); - - // cout << "AliGenUHKM::Generate() mother index in fParticles = " - // << (iparticle->GetFirstMother()==-1 ? -1 : iparticle->GetFirstMother()+1) - // << " ; hasMother = " << hasMother << endl; - // cout << "AliGenUHKM::Generate() type (0-hydro;1-jet) = " << iparticle->GetStatusCode() << endl; Bool_t hasDaughter = (iparticle->GetNDaughters() > 0); - - // cout << "AliGenUHKM::Generate() n.daughters = " << iparticle->GetNDaughters() - //<< " ; hasDaughter = " << hasDaughter << endl; - - + if(hasDaughter) { - // cout << "AliGenUHKM::Generate() decayed particle (not trackable)" << endl; // This particle has decayed // It will not be tracked // Add it only once with coordinates not // smeared with primary vertex position - Float_t p[3] = {p[0] = iparticle->Px(), - p[1] = iparticle->Py(), - p[2] = iparticle->Pz()}; + + // rotate the direction of the particle + partMomPhi = TMath::ATan2(iparticle->Py(), iparticle->Px()); + partPt = TMath::Hypot(iparticle->Px(), iparticle->Py()); + p[0] = partPt*TMath::Cos(partMomPhi+eventRotation); + p[1] = partPt*TMath::Sin(partMomPhi+eventRotation); + 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]); - v[0] = iparticle->Vx(); - v[1] = iparticle->Vy(); - v[2] = iparticle->Vz(); - Float_t time = iparticle->T(); - // if(kf==130 || kf==-130) cout << "type = " << iparticle->GetStatusCode() << endl; //j1/h0 + // rotate the freezeout point + partVtxPhi = TMath::ATan2(iparticle->Vy(), iparticle->Vx()); + partVtxR = TMath::Hypot(iparticle->Vx(), iparticle->Vy()); + v[0] = partVtxR*TMath::Cos(partVtxPhi + eventRotation); + v[1] = partVtxR*TMath::Cos(partVtxPhi + eventRotation); + v[2] = iparticle->Vz(); + Float_t time = iparticle->T(); + Int_t imo = -1; if(hasMother) { - imo = iparticle->GetFirstMother(); //index of mother particle in fParticles + imo = iparticle->GetFirstMother(); //index of mother particle in fParticles } // if has mother Bool_t trackFlag = kFALSE; // tFlag is kFALSE --> do not track the particle - - // printf("Pushing Track %d with status %d mother %d\n", kf, trackFlag, imo>=0?idsOnStack[imo]:imo); - PushTrack(trackFlag, (imo>=0 ? idsOnStack[imo+1] : imo), kf, - p[0], p[1], p[2], energy, - v[0], v[1], v[2], time, - polar[0], polar[1], polar[2], - (hasMother ? kPDecay : kPNoProcess), nt); - // cout << "AliGenUHKM::Generate() pushed on stack with stack index = " << nt - // << "; mother index on stack = " << (imo>=0 ? idsOnStack[imo+1] : imo) << endl; + + PushTrack(trackFlag, (imo>=0 ? idsOnStack[imo] : imo), kf, + p[0], p[1], p[2], energy, + v[0], v[1], v[2], time, + polar[0], polar[1], polar[2], + (hasMother ? kPDecay : kPNoProcess), nt); idsOnStack[i] = nt; + fNprimaries++; KeepTrack(nt); } else { - // cout << "AliGenUHKM::Generate() final particle --> push it twice on the stack" << endl; // This is a final state particle // It will be tracked // Add it TWICE to the stack !!! @@ -429,58 +411,64 @@ void AliGenUHKM::Generate() // this one will not be tracked // Second time with event-wide c0ordinates and vertex smearing // this one will be tracked - Float_t p[3] = {p[0] = iparticle->Px(), - p[1] = iparticle->Py(), - p[2] = iparticle->Pz()}; - mass = TDatabasePDG::Instance()->GetParticle(kf)->Mass(); + + // rotate the direction of the particle + partMomPhi = TMath::ATan2(iparticle->Py(), iparticle->Px()); + partPt = TMath::Hypot(iparticle->Px(), iparticle->Py()); + p[0] = partPt*TMath::Cos(partMomPhi+eventRotation); + p[1] = partPt*TMath::Sin(partMomPhi+eventRotation); + p[2] = iparticle->Pz(); + energy = sqrt(mass*mass + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]); - v[0] = iparticle->Vx(); - v[1] = iparticle->Vy(); - v[2] = iparticle->Vz(); + // rotate the freezeout point + partVtxPhi = TMath::ATan2(iparticle->Vy(), iparticle->Vx()); + partVtxR = TMath::Hypot(iparticle->Vx(), iparticle->Vy()); + v[0] = partVtxR*TMath::Cos(partVtxPhi + eventRotation); + v[1] = partVtxR*TMath::Cos(partVtxPhi + eventRotation); + v[2] = iparticle->Vz(); + Int_t type = iparticle->GetStatusCode(); // 1-from jet / 0-from hydro - //if(kf==130 || kf==-130) cout << "type = " << iparticle->GetStatusCode() << endl; //j1/h0 Int_t coeffT=1; if(type==1) coeffT=-1; //to separate particles from jets - + + Int_t imo = -1; if(hasMother) { - imo = iparticle->GetFirstMother(); + imo = iparticle->GetFirstMother(); } // if has mother + Bool_t trackFlag = kFALSE; // tFlag = kFALSE --> do not track this one, its for femtoscopy - PushTrack(trackFlag, (imo>=0 ? idsOnStack[imo+1] : imo), kf, - p[0], p[1], p[2], energy, - v[0], v[1], v[2], (iparticle->T())*coeffT, - polar[0], polar[1], polar[2], - hasMother ? kPDecay:kPNoProcess, nt); - // cout << "AliGenUHKM::Generate() pushed on stack with stack index = " << nt - // << "; mother index on stack = " << (imo>=0 ? idsOnStack[imo+1] : imo) << endl; - + PushTrack(trackFlag, (imo>=0 ? idsOnStack[imo] : imo), kf, + p[0], p[1], p[2], energy, + v[0], v[1], v[2], (iparticle->T())*coeffT, // freeze-out time is negative if the particle comes from jet + polar[0], polar[1], polar[2], + hasMother ? kPDecay:kPNoProcess, nt); + idsOnStack[i] = nt; fNprimaries++; KeepTrack(nt); - + origin[0] = origin0[0]+v[0]; origin[1] = origin0[1]+v[1]; origin[2] = origin0[2]+v[2]; + Float_t time = time0+iparticle->T(); imo = nt; - trackFlag = kTRUE; // tFlag = kTRUE --> track this one - // cout << "AliGenUHKM::Generate() trackFlag = " << trackFlag << endl; - + trackFlag = fTrackIt; // Track this particle, unless otherwise specified by fTrackIt + PushTrack(trackFlag, 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); - // cout << "AliGenUHKM::Generate() pushed on stack with stack index = " << nt - // << "; mother index on stack = " << imo << endl; + p[0], p[1], p[2], energy, + origin[0], origin[1], origin[2], time, + polar[0], polar[1], polar[2], + hasMother ? kPDecay:kPNoProcess, nt); + fNprimaries++; KeepTrack(nt); } } - + SetHighWaterMark(fNprimaries); TArrayF eventVertex; @@ -488,30 +476,36 @@ void AliGenUHKM::Generate() eventVertex[0] = origin0[0]; eventVertex[1] = origin0[1]; eventVertex[2] = origin0[2]; + Float_t eventTime = time0; // Builds the event header, to be called after each event AliGenEventHeader* header = new AliGenHijingEventHeader("UHKM"); + Double_t b = 0.; + Double_t npart = 0; + Double_t nbin = 0; + fUHKMgen->GetCentrality(b, npart, nbin); + printf("********** %13.3f %13.3f %13.3f \n", b, npart, nbin); + + ((AliGenHijingEventHeader*) header)->SetNProduced(fNprimaries); ((AliGenHijingEventHeader*) header)->SetPrimaryVertex(eventVertex); - ((AliGenHijingEventHeader*) header)->SetImpactParameter(0.0); + ((AliGenHijingEventHeader*) header)->SetInteractionTime(eventTime); + ((AliGenHijingEventHeader*) header)->SetImpactParameter(b); ((AliGenHijingEventHeader*) header)->SetTotalEnergy(0.0); ((AliGenHijingEventHeader*) header)->SetHardScatters(0); - ((AliGenHijingEventHeader*) header)->SetParticipants(0, 0); - ((AliGenHijingEventHeader*) header)->SetCollisions(0, 0, 0, 0); + ((AliGenHijingEventHeader*) header)->SetParticipants(Int_t(npart), 0); + ((AliGenHijingEventHeader*) header)->SetCollisions(Int_t(nbin), 0, 0, 0); ((AliGenHijingEventHeader*) header)->SetSpectators(0, 0, 0, 0); ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(0);//evrot); header->SetPrimaryVertex(fVertex); + header->SetInteractionTime(fTime); AddHeader(header); fCollisionGeometry = (AliGenHijingEventHeader*) header; - delete idsOnStack; - - // gAlice->SetGenEventHeader(header); - - // printf(" Finish Generate .. %d ..\n",nt); - // cout << "AliGenUHKM::Generate() OUT" << endl; + delete [] idsOnStack; + delete [] newPos; } void AliGenUHKM::Copy(TObject &) const @@ -522,8 +516,6 @@ void AliGenUHKM::Copy(TObject &) const void AliGenUHKM::SetAllParameters() { // Forward all input parameters to the TGenerator::TUHKMgen object - // cout << "AliGenUHKM::SetAllParameters() IN" << endl; - fUHKMgen->SetEcms(fHydjetParams.fSqrtS); fUHKMgen->SetBmin(fHydjetParams.fBmin); fUHKMgen->SetBmax(fHydjetParams.fBmax); @@ -546,7 +538,6 @@ void AliGenUHKM::SetAllParameters() { fUHKMgen->SetCoordAsymmPar(fHydjetParams.fEpsilon); fUHKMgen->SetGammaS(fHydjetParams.fCorrS); - // fUHKMgen->SetHdec(fHydjetParams.fTime); fUHKMgen->SetEtaType(fHydjetParams.fEtaType); fUHKMgen->SetFlagWeakDecay(fHydjetParams.fWeakDecay); @@ -569,17 +560,17 @@ void AliGenUHKM::SetAllParameters() { // fUHKMgen->SetMinimumMass(fMinMass); // fUHKMgen->SetMaximumMass(fMaxMass); - cout << "AliGenUHKM::Init() no. stable flagged particles = " << fStableFlagged << endl; + // cout << "AliGenUHKM::Init() no. stable flagged particles = " << fStableFlagged << endl; for(Int_t i=0; iPDGInfo(); // UHKM's PDG table TParticlePDG *rootTestParticle; ParticlePDG *uhkmTestParticle; - // cout << "particles with good status in UHKM table = " << uhkmPDG->GetNParticles() << endl; // loop over all particles in the SHARE table for(Int_t i=0; iGetNParticles(); i++) { - // cout << "particle #" << i << " ================" << endl; // get a particle specie uhkmTestParticle = uhkmPDG->GetPDGParticleByIndex(i); - // cout << "PDG = " << uhkmTestParticle->GetPDG() << endl; // check if this code exists in ROOT's table rootTestParticle = TDatabasePDG::Instance()->GetParticle(uhkmTestParticle->GetPDG()); if(!rootTestParticle) { // if not then add it to the ROOT's PDG database @@ -660,14 +645,10 @@ void AliGenUHKM::CheckPDGTable() { uhkmTestParticle->GetWidth(), (Int_t(uhkmTestParticle->GetBaryonNumber())==0 ? "meson" : "baryon"), uhkmTestParticle->GetPDG()); - // cout << "Not found in ROOT's PDG database --> added now" << endl; if(uhkmTestParticle->GetWidth()<1e-10) - cout << uhkmTestParticle->GetPDG() << "its mass " + cout << uhkmTestParticle->GetPDG() << " with mass " << TDatabasePDG::Instance()->GetParticle(uhkmTestParticle->GetPDG())->Mass() << TDatabasePDG::Instance()->GetParticle(uhkmTestParticle->GetPDG())->Width() << endl; } - // else - // cout << "Found in ROOT's PDG database --> not added" << endl; } // end for - // cout << "AliGenUHKM::CheckPDGTable() OUT" << endl; }