StRoot: StJetMaker/StJetReader.cxx Source File

00001 //M.L. Miller
00002 //MIT Software
00003 //6/04
00004 
00005 //std
00006 #include <map>
00007 #include <string>
00008 #include <algorithm>
00009 #include <iostream>
00010 
00011 
00012 //StEmc
00013 #include "StEmcClusterCollection.h"
00014 #include "StEmcPoint.h"
00015 #include "StEmcUtil/geometry/StEmcGeom.h"
00016 #include "StEmcUtil/others/emcDetectorName.h"
00017 #include "StEmcADCtoEMaker/StBemcData.h"
00018 #include "StEmcADCtoEMaker/StEmcADCtoEMaker.h"
00019 
00020 //root
00021 #include "TTree.h"
00022 #include "TFriendElement.h"
00023 #include "TFile.h"
00024 
00025 //StMuDst
00026 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00027 #include "StMuDSTMaker/COMMON/StMuDst.h"
00028 #include "StMuDSTMaker/COMMON/StMuEvent.h"
00029 
00030 //StJetMaker
00031 #include "StSpinPool/StJets/StJet.h"
00032 #include "StSpinPool/StJets/StJets.h"
00033 #include "StJetMaker/StJetReader.h"
00034 #include "StSpinPool/StJetSkimEvent/StJetSkimEvent.h"
00035 
00036 ClassImp(StJetReader);
00037 
00038 StJetReader::StJetReader(const char* name)
00039   : StMaker(name)
00040   , mFile(0)
00041   , mTree(0)
00042   , mSkimFile(0)
00043   , mSkimTree(0)
00044   , mValid(false)
00045   , mSkimEvent(0)
00046   , mOfstream(0)
00047 {
00048   LOG_INFO << "StJetReader::StJetReader()" << endm;
00049 }
00050 
00051 StJetReader::~StJetReader()
00052 {
00053   LOG_INFO << "StJetReader::~StJetReader()" << endm;
00054 }
00055 
00056 Int_t StJetReader::Init()
00057 {
00058   return kStOk;
00059 }
00060 
00061 void StJetReader::InitFile(const char* file)
00062 {
00063   LOG_INFO << "StJetReader::InitFile()" << endm;
00064 
00065   LOG_INFO << "open file:\t" << file << "\tfor reading" << endm;
00066   mFile = TFile::Open(file);
00067   assert(mFile);
00068 
00069   LOG_INFO << "recover tree" << endm;
00070   mTree = dynamic_cast<TTree*>(mFile->Get("jet"));
00071   assert(mTree);
00072 
00073   // Build index table to access jets by run number and event number
00074   int nentries = mTree->BuildIndex("mRunNumber","mEventNumber");
00075   if (nentries < 0) {
00076     LOG_WARN << "Could not build jet tree index: nentries=" << nentries << endm;
00077   }
00078   else {
00079     LOG_INFO << "Built jet tree index: nentries=" << nentries << endm;
00080   }
00081         
00082   LOG_INFO << "\tset tree pointer" << endm;
00083   LOG_INFO << "Number of entries in tree:\t" << mTree->GetEntriesFast() << endm;
00084 
00085   LOG_INFO << "\tGet Branches" << endm;
00086   TObjArray* branches = mTree->GetListOfBranches();
00087   if (!branches) {
00088     LOG_ERROR << "StJetReader::InitFile().  Null branches" << endm;
00089     abort();
00090   }
00091 
00092   LOG_INFO << "\tLoop on branches" << endm;
00093 
00094   for (int i=0; i<branches->GetEntriesFast(); ++i) {
00095     TBranch* branch = dynamic_cast<TBranch*>(branches->At(i));
00096     if (!branch) {
00097       LOG_ERROR << "StJetReader::InitFile().  Null branch" << endm;
00098       abort();
00099     }
00100 
00101     LOG_INFO << "\t--- Found branch:\t" << branch->GetName() << endm;
00102   }
00103 
00104   if (0) {
00105     string jetCheck(file);
00106     jetCheck += ".read.txt";
00107     mOfstream = new ofstream(jetCheck.c_str());
00108   }
00109 
00110   LOG_INFO << "\tfinished!" <<endm;
00111 
00112   return;
00113 }
00114 
00115 void StJetReader::InitJetSkimFile(const char* file)
00116 {
00117   LOG_INFO << "StJetReader::InitJetSkimFile(const char* file)" << endm;
00118   LOG_INFO << "open file:\t" << file << "\tfor reading" << endm;
00119 
00120   mSkimFile = TFile::Open(file);
00121   assert(mSkimFile);
00122 
00123   mSkimTree = dynamic_cast<TTree*>(mSkimFile->Get("jetSkimTree"));
00124   assert(mSkimTree);
00125 
00126   // Build index table to access skim events by run number and event number
00127   int nentries = mSkimTree->BuildIndex("mRunId","mEventId");
00128   if (nentries < 0) {
00129     LOG_WARN << "Could not build skim tree index: nentries=" << nentries << endm;
00130   }
00131   else {
00132     LOG_INFO << "Built skim tree index: nentries=" << nentries << endm;
00133   }
00134 
00135   mSkimEvent = new StJetSkimEvent;
00136 
00137   //set branch status
00138   mSkimTree->SetBranchStatus("skimEventBranch", 1);
00139   mSkimTree->SetBranchAddress("skimEventBranch", &mSkimEvent);
00140 
00141   LOG_INFO << "successfully recovered tree from file." << endm;
00142 }
00143 
00144 void StJetReader::InitTree(TTree* tree)
00145 {
00146   LOG_INFO << "StJetReader::InitTree()" << endm;
00147 
00148   LOG_INFO << "\tset tree pointer" << endm;
00149   LOG_INFO << "Number of entries in tree:\t" << tree->GetEntriesFast() << endm;
00150 
00151   TList* friendList = tree->GetListOfFriends();
00152   TTree* t=0;
00153   for (int j=0; j<friendList->GetSize(); ++j) {
00154     TFriendElement* fr = static_cast<TFriendElement*>( friendList->At(j) );
00155     string tree_name( fr->GetTreeName() );
00156     if (tree_name == "jet") {
00157       t = fr->GetTree();
00158       break;
00159     }
00160   }
00161   assert(t);
00162 
00163   LOG_INFO << "\tGet Branches" << endm;
00164   TObjArray* branches = t->GetListOfBranches();
00165   if (!branches) {
00166     LOG_ERROR << "StJetReader::InitFile().  Null branches" << endm;
00167     abort();
00168   }
00169 
00170   LOG_INFO << "\tLoop on branches" << endm;
00171 
00172   for (int i=0; i<branches->GetEntriesFast(); ++i) {
00173     TBranch* branch = dynamic_cast<TBranch*>(branches->At(i));
00174     if (!branch) {
00175       LOG_INFO << "StJetReader::InitFile().  Null branch" << endm;
00176       abort();
00177     }
00178 
00179     LOG_INFO << "\t--- Found branch:\t" << branch->GetName() << endm;
00180   }
00181 
00182   LOG_INFO << "\tfinished!" << endm;
00183 
00184   return;
00185 }
00186 
00187 int StJetReader::preparedForDualRead()
00188 {
00189   LOG_INFO << "StJetReader::preparedForDualRead()" << endm;
00190   int jetStatus = mTree != 0;
00191   int skimStatus = mSkimTree != 0;
00192   LOG_INFO << "jet tree status:\t" << jetStatus << endm;
00193   LOG_INFO << "skim tree status:\t" << skimStatus << endm;
00194   assert(mTree);
00195   assert(mSkimTree);
00196 
00197   int nJetEvents = mTree->GetEntriesFast();
00198   int nSkimEvents = mSkimTree->GetEntriesFast();
00199   LOG_INFO << "nJetEvents:\t" << nJetEvents << endm;
00200   LOG_INFO << "nSkimEvents:\t" << nSkimEvents << endm;
00201   assert(nJetEvents == nSkimEvents);
00202   assert(mSkimEvent);
00203 
00204   LOG_INFO << "Congratulations, you are ready to process events!" << endm;
00205   mValid = true;
00206   return 1;
00207 }
00208 
00209 Int_t StJetReader::Make()
00210 {
00211   if (mTree && mTree->GetEntryWithIndex(GetRunNumber(),GetEventNumber()) < 0) {
00212     LOG_WARN << Form("Could not find jet tree entry for run=%d event=%d",GetRunNumber(),GetEventNumber()) << endm;
00213     return kStWarn;
00214   }
00215 
00216   if (mSkimTree && mSkimTree->GetEntryWithIndex(GetRunNumber(),GetEventNumber()) < 0) {
00217     LOG_WARN << Form("Could not find skim tree entry for run=%d event=%d",GetRunNumber(),GetEventNumber()) << endm;
00218     return kStWarn;
00219   }
00220 
00221   return kStOk;
00222 }
00223 
00224 Int_t StJetReader::Finish()
00225 {
00226   if (mOfstream) {
00227     delete mOfstream;
00228     mOfstream=0;
00229   }
00230         
00231   return kStOk;
00232 }
00233 
00234 void dumpProtojetToStream(int event, ostream& os, StJets* stjets);
00235 
00236 string idString(TrackToJetIndex* t)
00237 {
00238   string idstring;
00239   StDetectorId mDetId = t->detectorId();
00240   if (mDetId==kTpcId) {
00241     idstring = "kTpcId";
00242   }
00243   else if (mDetId==kBarrelEmcTowerId) {
00244     idstring = "kBarrelEmcTowerId";
00245   }
00246   else if (mDetId==kEndcapEmcTowerId) {
00247     idstring = "kEndcapEmcTowerId";
00248   }
00249   else {
00250     idstring = "kUnknown";
00251   }
00252   return idstring;
00253 }
00254 
00255 //nice check to verify that the jet 4-mom is equal to the _vector_ sum of it's part
00256 bool verifyJet(StJets* stjets, int ijet)
00257 {
00258   TClonesArray& jets = *(stjets->jets());
00259   StJet* pj = static_cast<StJet*>(jets[ijet]);
00260     
00261     
00262   StThreeVectorD j3(pj->Px(), pj->Py(), pj->Pz());
00263   StLorentzVectorD j4(pj->E(),j3 );
00264         
00265   //LOG_INFO <<"\tjet:\t"<<ijet<<"\t4mom:\t"<<j4<<endm;
00266         
00267   StLorentzVectorD jetMom(0., 0., 0., 0.);
00268     
00269   typedef vector<TLorentzVector*> FourpList;
00270   FourpList particles = stjets->particles(ijet);
00271     
00272   for (FourpList::iterator it=particles.begin(); it!=particles.end(); ++it) {
00273     TLorentzVector* v = *it;
00274     StThreeVectorD v3(v->Px(), v->Py(), v->Pz() );
00275     StLorentzVectorD v4(v->E(), v3 );
00276     jetMom += v4;
00277     //LOG_INFO <<"\t\t4p:\t"<<v4<<endm;
00278   }
00279   //LOG_INFO <<"\t\t\tcheck:\t"<<jetMom<<endm;
00280   StLorentzVectorD diff = j4-jetMom;
00281   if (abs(diff)>1.e-6) { //they have to be the same to 1 eV
00282     LOG_INFO << "verifyJet.  assert will fail for jet:\t" << ijet << "\t4p:\t" << j4 <<"\tcompared to sum_particles:\t" << jetMom << endm;
00283     return false;
00284   }
00285   else {
00286     return true;
00287   }
00288 }
00289 
00290 
00291 void StJetReader::exampleFastAna()
00292 {
00293   LOG_INFO << "StJetReader::exampleEventAna()" << endm;
00294         
00295   StJetSkimEvent* skEv = skimEvent();
00296             
00297   //basic information:
00298   LOG_INFO << "fill/run/event:\t" << skEv->fill() << "\t" << skEv->runId() << "\t" << skEv->eventId() << endm;
00299   LOG_INFO << "fileName:\t" << skEv->mudstFileName().GetString() << endm;
00300         
00301   //some spin info:
00302   LOG_INFO << "bx48:\t" << skEv->bx48() << "\tspinBits:\t" << skEv->spinBits() << endm;
00303         
00304   //some event info:
00305   LOG_INFO << "Ebbc:\t" << skEv->eBbc() << "\tWbbc:\t" << skEv->wBbc() << "\tbbcTimeBin:\t" << skEv->bbcTimeBin() << endm;
00306         
00307   //best vertex info:
00308   const float* bestPos = skEv->bestVert()->position();
00309   LOG_INFO << "best Vertex (x,y,z):\t" << bestPos[0] << "\t" << bestPos[1] << "\t" << bestPos[2] << endm;
00310         
00311   //trigger info:
00312   const TClonesArray* trigs = skEv->triggers();
00313   assert(trigs);
00314   int nTrigs = trigs->GetEntriesFast();
00315   for (int i=0; i<nTrigs; ++i) {
00316     StJetSkimTrig* trig = static_cast<StJetSkimTrig*>( trigs->At(i) );
00317     assert(trig);
00318     if (trig->didFire() != 0) {
00319       //LOG_INFO << "\tTriggerd:\t" << trig->trigId() << "\twith prescale:\t" << trig->prescale() << endm;
00320     }
00321     if (trig->shouldFireL2() == 1) {
00322       LOG_INFO << "\tshTrigger L2: " << trig->trigId() << endm;
00323       const int* l2temp = trig->L2ResultEmulated();
00324       for (int ii = 0; ii < 9; ++ii) {
00325         LOG_INFO << "SimL2--- " << ii << "\t" << l2temp[ii] << endm;
00326       }
00327     }
00328   }
00329                 
00330   //L2 Info:
00331   LOG_INFO << "\n--- Non-zero L2 Results:" << endm;
00332   const int* l2temp = skEv->L2Result();
00333   for (int ii=0; ii<9; ii++) {
00334     if (l2temp[ii]!=0)  LOG_INFO << "DatL2--- " << ii << "\t" << l2temp[ii] << endm;
00335   }
00336 
00337   //vertex info:
00338   const TClonesArray* verts = skEv->vertices();
00339   assert(verts);
00340   int nVerts = verts->GetEntriesFast();
00341   for (int i=0; i<nVerts; ++i) {
00342     StJetSkimVert* vert = static_cast<StJetSkimVert*>( verts->At(i) );
00343     assert(vert);
00344     const float* vertPos = vert->position();
00345     LOG_INFO << "vert:\t" << i << "\t at z:\t" << vertPos[2] << "\tranking:\t"
00346              << vert->ranking() << "\tnTracks:\t" << vert->nTracksUsed() << endm;
00347   }
00348 
00349   //And then access to jets from different algorithms...
00350   LOG_INFO << "\nLoop on jet branches" << endm;
00351   for (int i = 0; i < numberOfBranches(); ++i) {
00352     StJets* stjets = getStJets(i);
00353     int nJets = stjets->nJets();
00354     LOG_INFO << "Found\t" << nJets << "\tjets from:\t" << branchName(i) << endm;
00355 
00356     //first make sure that we have the same run/event:
00357     assert(stjets->runId()==skEv->runId());
00358     assert(stjets->eventId()==skEv->eventId());
00359 
00360     TClonesArray* jets = stjets->jets();
00361                 
00362     for(int ijet=0; ijet<nJets; ++ijet){ 
00363 
00364       //loop on jets
00365       StJet* j = dynamic_cast<StJet*>( jets->At(ijet) );
00366       assert(j);
00367       //assert(verifyJet(stjets, ijet));
00368                         
00369       LOG_INFO <<"jet:\t"<<ijet<<"\tEjet:\t"<<j->E()<<"\tEta:\t"<<j->Eta()<<"\tPhi:\t"<<j->Phi()<<"\tdetEta:\t"<<j->detEta()<<endm;
00370                         
00371       //look at 4-momenta in the jet:
00372       typedef vector<TLorentzVector*> TrackToJetVec;
00373       TrackToJetVec particles = stjets->particles(ijet);
00374                         
00375       for (TrackToJetVec::iterator partIt=particles.begin(); partIt!=particles.end(); ++partIt) {
00376         const TLorentzVector* theParticle = *partIt;
00377         LOG_INFO <<"\tparticle \t pt:\t"<<theParticle->Pt()<<"\teta:\t"<<theParticle->Eta()<<"\tphi:\t"<<theParticle->Phi()<<endm;
00378       }
00379     }
00380   }     
00381 }
00382 
00383 void StJetReader::exampleEventAna()
00384 {
00385   LOG_INFO <<"StJetReader::exampleEventAna()"<<endm;
00386         
00387   LOG_INFO <<"nPrimary:\t"<<StMuDst::numberOfPrimaryTracks() << endm;
00388 
00389   for (int i = 0; i < numberOfBranches(); ++i) {
00390     StJets* stjets = getStJets(i);
00391     int nJets = stjets->nJets();
00392     LOG_INFO <<"Found\t"<<nJets<<"\tjets from:\t"<<branchName(i)<<endm;
00393                 
00394     TClonesArray* jets = stjets->jets();
00395                 
00396     //Dylan, here's a nice check...
00397     if (0) {
00398       if (stjets->nJets()>0) {
00399         dumpProtojetToStream(StMuDst::event()->eventId(), *mOfstream, stjets);
00400       }
00401     }
00402                 
00403     for(int ijet=0; ijet<nJets; ++ijet){ 
00404                         
00405       //loop on jets
00406       StJet* j = dynamic_cast<StJet*>( jets->At(ijet) );
00407       assert(j);
00408       assert(verifyJet(stjets, ijet));
00409                         
00410       LOG_INFO <<"jet:\t"<<ijet<<"\tEjet:\t"<<j->E()<<"\tEta:\t"<<j->Eta()<<"\tPhi:\t"<<j->Phi()<<endm;
00411                         
00412       //look at 4-momenta in the jet:
00413       typedef vector<TLorentzVector*> TrackToJetVec;
00414       TrackToJetVec particles = stjets->particles(ijet);
00415     }
00416   }
00417 }
00418 
00419 int StJetReader::numberOfBranches() const
00420 {
00421   return mTree->GetNbranches();
00422 }
00423 
00424 const char* StJetReader::branchName(int i) const
00425 {
00426   TBranch* branch = (TBranch*)mTree->GetListOfBranches()->At(i);
00427   return branch ? branch->GetName() : 0;
00428 }
00429 
00430 StJets* StJetReader::getStJets(int i) const
00431 {
00432   TBranch* branch = (TBranch*)mTree->GetListOfBranches()->At(i);
00433   return branch ? *(StJets**)branch->GetAddress() : 0;
00434 }
00435 
00436 StJets* StJetReader::getStJets(const char* bname) const
00437 {
00438   TBranch* branch = mTree->GetBranch(bname);
00439   return branch ? *(StJets**)branch->GetAddress() : 0;
00440 }