StRoot: StSpinPool/StBfcTriggerFilterMaker/StPythiaEventMaker.cxx Source File

00001 //
00002 // Pibero Djawotho <pibero@tamu.edu>
00003 // Texas A&M University
00004 // 4 May 2010
00005 //
00006 
00007 // ROOT
00008 #include "TFile.h"
00009 #include "TTree.h"
00010 
00011 // STAR
00012 #include "StSpinPool/StJetSkimEvent/StPythiaEvent.h"
00013 #include "tables/St_g2t_event_Table.h"
00014 #include "tables/St_g2t_pythia_Table.h"
00015 #include "tables/St_g2t_vertex_Table.h"
00016 #include "tables/St_particle_Table.h"
00017 #include "StSpinPool/StMCAsymMaker/StMCAsymMaker.h"
00018 
00019 // Local
00020 #include "StPythiaEventMaker.h"
00021 
00022 ClassImp(StPythiaEventMaker);
00023 
00024 void StPythiaEventMaker::Clear(Option_t* option)
00025 {
00026   mPythiaEvent->Clear(option);
00027 }
00028 
00029 int StPythiaEventMaker::Init()
00030 {
00031   assert(!mFileName.IsNull());
00032   mFile = TFile::Open(mFileName,"recreate");
00033   assert(mFile);
00034   mPythiaEvent = new StPythiaEvent;
00035   mTree = new TTree("PythiaTree","Pythia Record");
00036   mTree->Branch("PythiaBranch","StPythiaEvent",&mPythiaEvent);
00037   return kStOk;
00038 }
00039 
00040 int StPythiaEventMaker::Make()
00041 {
00042   getEvent();
00043   getPythia();
00044   getVertex();
00045   getParticles();
00046   getAsymmetries();
00047 
00048   if (Debug()) { LOG_DEBUG << *mPythiaEvent << endm; }
00049 
00050   mTree->Fill();
00051 
00052   return kStOk;
00053 }
00054 
00055 int StPythiaEventMaker::Finish()
00056 {
00057   mFile->Write();
00058   mFile->Close();
00059   return kStOk;
00060 }
00061 
00062 void StPythiaEventMaker::getEvent()
00063 {
00064 //   TDataSet* geant = GetDataSet("geant");
00065 //   if (geant) {
00066 //     TDataSetIter iter(geant);
00067 //     St_g2t_event* eventDescriptor = (St_g2t_event*)iter("g2t_event");
00068 //     if (eventDescriptor) {
00069 //       g2t_event_st* eventTable = (g2t_event_st*)eventDescriptor->GetTable();
00070 //       if (eventTable) {
00071 //      mPythiaEvent->setRunId(eventTable->n_run);
00072 //      mPythiaEvent->setEventId(eventTable->n_event);
00073 //       }
00074 //     }
00075 //   }
00076 
00077   mPythiaEvent->setRunId(GetRunNumber());
00078   mPythiaEvent->setEventId(GetEventNumber());
00079 }
00080 
00081 void StPythiaEventMaker::getPythia()
00082 {
00083   TDataSet* geant = GetDataSet("geant");
00084   if (geant) {
00085     TDataSetIter iter(geant);
00086     St_g2t_pythia* pythiaDescriptor = (St_g2t_pythia*)iter("g2t_pythia");
00087     if (pythiaDescriptor) {
00088       g2t_pythia_st* pythiaTable = (g2t_pythia_st*)pythiaDescriptor->GetTable();
00089       if (pythiaTable) {
00090         mPythiaEvent->setProcessId(pythiaTable->subprocess_id);
00091         mPythiaEvent->setS(pythiaTable->mand_s);
00092         mPythiaEvent->setT(pythiaTable->mand_t);
00093         mPythiaEvent->setU(pythiaTable->mand_u);
00094         mPythiaEvent->setPt(pythiaTable->hard_p);
00095         mPythiaEvent->setCosTheta(pythiaTable->cos_th);
00096         mPythiaEvent->setX1(pythiaTable->bjor_1);
00097         mPythiaEvent->setX2(pythiaTable->bjor_2);
00098         mPythiaEvent->setMstu72(pythiaTable->mstu72);
00099         mPythiaEvent->setMstu73(pythiaTable->mstu73);
00100         mPythiaEvent->setMstp111(pythiaTable->mstp111);
00101       }
00102     }
00103   }
00104 }
00105 
00106 void StPythiaEventMaker::getVertex()
00107 {
00108   TDataSet* geant = GetDataSet("geant");
00109   if (geant) {
00110     TDataSetIter iter(geant);
00111     St_g2t_vertex* vertexDescriptor = (St_g2t_vertex*)iter("g2t_vertex");
00112     if (vertexDescriptor) {
00113       g2t_vertex_st* vertexTable = (g2t_vertex_st*)vertexDescriptor->GetTable();
00114       if (vertexTable) {
00115         mPythiaEvent->setVertex(TVector3(vertexTable[0].ge_x));
00116       }
00117     }
00118   }
00119 }
00120 
00121 void StPythiaEventMaker::getParticles()
00122 {
00123   TDataSet* geant = GetDataSet("geant");
00124   if (geant) {
00125     TDataSetIter iter(geant);
00126     St_particle* particleDescriptor = (St_particle*)iter("particle");
00127     if (particleDescriptor) {
00128       particle_st* particleTable = (particle_st*)particleDescriptor->GetTable();
00129       if (particleTable) {
00130         for (int i = 0; i < particleDescriptor->GetNRows(); ++i) {
00131           mPythiaEvent->addParticle(TParticle(particleTable[i].idhep, // pdg
00132                                               particleTable[i].isthep, // status
00133                                               particleTable[i].jmohep[0], // mother1
00134                                               particleTable[i].jmohep[1], // mother2
00135                                               particleTable[i].jdahep[0], // daughter1
00136                                               particleTable[i].jdahep[1], // daughter2
00137                                               TLorentzVector(particleTable[i].phep), // momentum and energy
00138                                               TLorentzVector(particleTable[i].vhep))); // production vertex and time
00139         }
00140       }
00141     }
00142   }
00143 }
00144 
00145 void StPythiaEventMaker::getAsymmetries()
00146 {
00147   float s = mPythiaEvent->s();
00148   float t = mPythiaEvent->t();
00149   float u = mPythiaEvent->u();
00150   int pid = mPythiaEvent->processId();
00151   int flavor1 = mPythiaEvent->particle(4)->GetPdgCode();
00152   int flavor2 = mPythiaEvent->particle(5)->GetPdgCode();
00153   int flavor3 = mPythiaEvent->particle(6)->GetPdgCode();
00154   int flavor4 = mPythiaEvent->particle(7)->GetPdgCode();
00155   float x1 = mPythiaEvent->x1();
00156   float x2 = mPythiaEvent->x2();
00157   float Q2 = mPythiaEvent->Q2();
00158 
00159   mPythiaEvent->setDF1(StPythiaEvent::LO,StMCAsymMaker::get_polPDF_LO(flavor1,x1,Q2));
00160   mPythiaEvent->setDF1(StPythiaEvent::NLO,StMCAsymMaker::get_polPDF_NLO(flavor1,x1,Q2));
00161   mPythiaEvent->setDF1(StPythiaEvent::ZERO,StMCAsymMaker::get_polPDF_NLO_g0(flavor1,x1,Q2));
00162   mPythiaEvent->setDF1(StPythiaEvent::MAX,StMCAsymMaker::get_polPDF_NLO_gmax(flavor1,x1,Q2));
00163   mPythiaEvent->setDF1(StPythiaEvent::MIN,StMCAsymMaker::get_polPDF_NLO_gmin(flavor1,x1,Q2));
00164   mPythiaEvent->setDF1(StPythiaEvent::M015,StMCAsymMaker::get_polPDF_NLO_m015(flavor1,x1,Q2));
00165   mPythiaEvent->setDF1(StPythiaEvent::M030,StMCAsymMaker::get_polPDF_NLO_m030(flavor1,x1,Q2));
00166   mPythiaEvent->setDF1(StPythiaEvent::M045,StMCAsymMaker::get_polPDF_NLO_m045(flavor1,x1,Q2));
00167   mPythiaEvent->setDF1(StPythiaEvent::M060,StMCAsymMaker::get_polPDF_NLO_m060(flavor1,x1,Q2));
00168   mPythiaEvent->setDF1(StPythiaEvent::M075,StMCAsymMaker::get_polPDF_NLO_m075(flavor1,x1,Q2));
00169   mPythiaEvent->setDF1(StPythiaEvent::M090,StMCAsymMaker::get_polPDF_NLO_m090(flavor1,x1,Q2));
00170   mPythiaEvent->setDF1(StPythiaEvent::M105,StMCAsymMaker::get_polPDF_NLO_m105(flavor1,x1,Q2));
00171   mPythiaEvent->setDF1(StPythiaEvent::P030,StMCAsymMaker::get_polPDF_NLO_p030(flavor1,x1,Q2));
00172   mPythiaEvent->setDF1(StPythiaEvent::P045,StMCAsymMaker::get_polPDF_NLO_p045(flavor1,x1,Q2));
00173   mPythiaEvent->setDF1(StPythiaEvent::P060,StMCAsymMaker::get_polPDF_NLO_p060(flavor1,x1,Q2));
00174   mPythiaEvent->setDF1(StPythiaEvent::P070,StMCAsymMaker::get_polPDF_NLO_p070(flavor1,x1,Q2));
00175   mPythiaEvent->setDF1(StPythiaEvent::GS_NLOA,StMCAsymMaker::get_polPDF_NLO_GSA(flavor1,x1,Q2));
00176   mPythiaEvent->setDF1(StPythiaEvent::GS_NLOB,StMCAsymMaker::get_polPDF_NLO_GSB(flavor1,x1,Q2));
00177   mPythiaEvent->setDF1(StPythiaEvent::GS_NLOC,StMCAsymMaker::get_polPDF_NLO_GSC(flavor1,x1,Q2));
00178   mPythiaEvent->setDF1(StPythiaEvent::DSSV,StMCAsymMaker::get_polPDF_NLO_DSSV(flavor1,x1,Q2));
00179   mPythiaEvent->setDF1(StPythiaEvent::LSS1,StMCAsymMaker::get_polPDF_NLO_LSS1(flavor1,x1,Q2));
00180   mPythiaEvent->setDF1(StPythiaEvent::LSS2,StMCAsymMaker::get_polPDF_NLO_LSS2(flavor1,x1,Q2));
00181   mPythiaEvent->setDF1(StPythiaEvent::LSS3,StMCAsymMaker::get_polPDF_NLO_LSS3(flavor1,x1,Q2));
00182   mPythiaEvent->setDF1(StPythiaEvent::AAC1,StMCAsymMaker::get_polPDF_NLO_AAC1(flavor1,x1,Q2));
00183   mPythiaEvent->setDF1(StPythiaEvent::AAC2,StMCAsymMaker::get_polPDF_NLO_AAC2(flavor1,x1,Q2));
00184   mPythiaEvent->setDF1(StPythiaEvent::AAC3,StMCAsymMaker::get_polPDF_NLO_AAC3(flavor1,x1,Q2));
00185   mPythiaEvent->setDF1(StPythiaEvent::BB1,StMCAsymMaker::get_polPDF_NLO_BB1(flavor1,x1,Q2));
00186   mPythiaEvent->setDF1(StPythiaEvent::BB2,StMCAsymMaker::get_polPDF_NLO_BB2(flavor1,x1,Q2));
00187   mPythiaEvent->setDF1(StPythiaEvent::DNS1,StMCAsymMaker::get_polPDF_NLO_DNS1(flavor1,x1,Q2));
00188   mPythiaEvent->setDF1(StPythiaEvent::DNS2,StMCAsymMaker::get_polPDF_NLO_DNS2(flavor1,x1,Q2));
00189 
00190   mPythiaEvent->setDF2(StPythiaEvent::LO,StMCAsymMaker::get_polPDF_LO(flavor2,x2,Q2));
00191   mPythiaEvent->setDF2(StPythiaEvent::NLO,StMCAsymMaker::get_polPDF_NLO(flavor2,x2,Q2));
00192   mPythiaEvent->setDF2(StPythiaEvent::ZERO,StMCAsymMaker::get_polPDF_NLO_g0(flavor2,x2,Q2));
00193   mPythiaEvent->setDF2(StPythiaEvent::MAX,StMCAsymMaker::get_polPDF_NLO_gmax(flavor2,x2,Q2));
00194   mPythiaEvent->setDF2(StPythiaEvent::MIN,StMCAsymMaker::get_polPDF_NLO_gmin(flavor2,x2,Q2));
00195   mPythiaEvent->setDF2(StPythiaEvent::M015,StMCAsymMaker::get_polPDF_NLO_m015(flavor2,x2,Q2));
00196   mPythiaEvent->setDF2(StPythiaEvent::M030,StMCAsymMaker::get_polPDF_NLO_m030(flavor2,x2,Q2));
00197   mPythiaEvent->setDF2(StPythiaEvent::M045,StMCAsymMaker::get_polPDF_NLO_m045(flavor2,x2,Q2));
00198   mPythiaEvent->setDF2(StPythiaEvent::M060,StMCAsymMaker::get_polPDF_NLO_m060(flavor2,x2,Q2));
00199   mPythiaEvent->setDF2(StPythiaEvent::M075,StMCAsymMaker::get_polPDF_NLO_m075(flavor2,x2,Q2));
00200   mPythiaEvent->setDF2(StPythiaEvent::M090,StMCAsymMaker::get_polPDF_NLO_m090(flavor2,x2,Q2));
00201   mPythiaEvent->setDF2(StPythiaEvent::M105,StMCAsymMaker::get_polPDF_NLO_m105(flavor2,x2,Q2));
00202   mPythiaEvent->setDF2(StPythiaEvent::P030,StMCAsymMaker::get_polPDF_NLO_p030(flavor2,x2,Q2));
00203   mPythiaEvent->setDF2(StPythiaEvent::P045,StMCAsymMaker::get_polPDF_NLO_p045(flavor2,x2,Q2));
00204   mPythiaEvent->setDF2(StPythiaEvent::P060,StMCAsymMaker::get_polPDF_NLO_p060(flavor2,x2,Q2));
00205   mPythiaEvent->setDF2(StPythiaEvent::P070,StMCAsymMaker::get_polPDF_NLO_p070(flavor2,x2,Q2));
00206   mPythiaEvent->setDF2(StPythiaEvent::GS_NLOA,StMCAsymMaker::get_polPDF_NLO_GSA(flavor2,x2,Q2));
00207   mPythiaEvent->setDF2(StPythiaEvent::GS_NLOB,StMCAsymMaker::get_polPDF_NLO_GSB(flavor2,x2,Q2));
00208   mPythiaEvent->setDF2(StPythiaEvent::GS_NLOC,StMCAsymMaker::get_polPDF_NLO_GSC(flavor2,x2,Q2));
00209   mPythiaEvent->setDF2(StPythiaEvent::DSSV,StMCAsymMaker::get_polPDF_NLO_DSSV(flavor2,x2,Q2));
00210   mPythiaEvent->setDF2(StPythiaEvent::LSS1,StMCAsymMaker::get_polPDF_NLO_LSS1(flavor2,x2,Q2));
00211   mPythiaEvent->setDF2(StPythiaEvent::LSS2,StMCAsymMaker::get_polPDF_NLO_LSS2(flavor2,x2,Q2));
00212   mPythiaEvent->setDF2(StPythiaEvent::LSS3,StMCAsymMaker::get_polPDF_NLO_LSS3(flavor2,x2,Q2));
00213   mPythiaEvent->setDF2(StPythiaEvent::AAC1,StMCAsymMaker::get_polPDF_NLO_AAC1(flavor2,x2,Q2));
00214   mPythiaEvent->setDF2(StPythiaEvent::AAC2,StMCAsymMaker::get_polPDF_NLO_AAC2(flavor2,x2,Q2));
00215   mPythiaEvent->setDF2(StPythiaEvent::AAC3,StMCAsymMaker::get_polPDF_NLO_AAC3(flavor2,x2,Q2));
00216   mPythiaEvent->setDF2(StPythiaEvent::BB1,StMCAsymMaker::get_polPDF_NLO_BB1(flavor2,x2,Q2));
00217   mPythiaEvent->setDF2(StPythiaEvent::BB2,StMCAsymMaker::get_polPDF_NLO_BB2(flavor2,x2,Q2));
00218   mPythiaEvent->setDF2(StPythiaEvent::DNS1,StMCAsymMaker::get_polPDF_NLO_DNS1(flavor2,x2,Q2));
00219   mPythiaEvent->setDF2(StPythiaEvent::DNS2,StMCAsymMaker::get_polPDF_NLO_DNS2(flavor2,x2,Q2));
00220 
00221   mPythiaEvent->setF1(StPythiaEvent::LO,StMCAsymMaker::get_unpolPDF_LO(flavor1,x1,Q2));
00222   mPythiaEvent->setF1(StPythiaEvent::NLO,StMCAsymMaker::get_unpolPDF_NLO(flavor1,x1,Q2));
00223 
00224   mPythiaEvent->setF2(StPythiaEvent::LO,StMCAsymMaker::get_unpolPDF_LO(flavor2,x2,Q2));
00225   mPythiaEvent->setF2(StPythiaEvent::NLO,StMCAsymMaker::get_unpolPDF_NLO(flavor2,x2,Q2));
00226 
00227   mPythiaEvent->setPartonALL(StMCAsymMaker::getPartonicALL(s,t,u,pid,flavor1,flavor2,flavor3,flavor4));
00228 }
00229 
00230 ostream& operator<<(ostream& out, const TVector3& v)
00231 {
00232   return out << v.x() << '\t' << v.y() << '\t' << v.z();
00233 }
00234 
00235 ostream& operator<<(ostream& out, const TParticle& particle)
00236 {
00237   out << " name=" << particle.GetName()
00238       << " pdg=" << particle.GetPdgCode()
00239       << " sta=" << particle.GetStatusCode()
00240       << " mo1=" << particle.GetMother(0)
00241       << " mo2=" << particle.GetMother(1)
00242       << " da1=" << particle.GetDaughter(0)
00243       << " da2=" << particle.GetDaughter(1)
00244       << " pt=" << particle.Pt()
00245       << " pz=" << particle.Pz()
00246       << " eta=" << particle.Eta()
00247       << " phi=" << particle.Phi()
00248       << " e=" << particle.Energy();
00249   return out;
00250 }
00251 
00252 ostream& operator<<(ostream& out, const StPythiaEvent& pythiaEvent)
00253 {
00254   out << "Pythia Record\n"
00255       << "run:\t"      << pythiaEvent.runId()      << '\n'
00256       << "event:\t"    << pythiaEvent.eventId()    << '\n'
00257       << "vertex:\t"   << pythiaEvent.vertex()     << '\n'
00258       << "pid:\t"      << pythiaEvent.processId()  << '\n'
00259       << "s:\t"        << pythiaEvent.s()          << '\n'
00260       << "t:\t"        << pythiaEvent.t()          << '\n'
00261       << "u:\t"        << pythiaEvent.u()          << '\n'
00262       << "pt:\t"       << pythiaEvent.pt()         << '\n'
00263       << "cosTheta:\t" << pythiaEvent.cosTheta()   << '\n'
00264       << "x1:\t"       << pythiaEvent.x1()         << '\n'
00265       << "x2:\t"       << pythiaEvent.x2()         << '\n';
00266   out << Form("%5s%14s%7s%9s%7s%9s%9s%9s%9s%9s\n",
00267               "line","particle/jet","status","PDG","mother","pt","eta","phi","E","m");
00268   for (int i = 0; i < pythiaEvent.numberOfParticles(); ++i) {
00269     const TParticle* particle = pythiaEvent.particle(i);
00270     out << Form("%5d%14s%7d%9d%7d%9.3f%9.3f%9.3f%9.3f%9.3f\n",
00271                 i+1,particle->GetName(),particle->GetStatusCode(),particle->GetPdgCode(),particle->GetFirstMother(),
00272                 particle->Pt(),particle->Eta(),particle->Phi(),particle->Energy(),particle->GetCalcMass());
00273   }
00274   return out;
00275 }