00001
00002
00003
00004
00005
00006 #include <map>
00007 #include <string>
00008 #include <algorithm>
00009 #include <iostream>
00010 #include <math.h>
00011
00012
00013 #include "TTree.h"
00014 #include "TFriendElement.h"
00015 #include "TFile.h"
00016
00017
00018 #include "StMessMgr.h"
00019
00020
00021 #include "StSpinPool/StMCAsymMaker/StMCAsymMaker.h"
00022 #include "StSpinPool/StMCAsymMaker/StPythiaEvent.h"
00023
00024
00025 #include "StMuDSTMaker/COMMON/StMuDstMaker.h"
00026
00027
00028 #include "StJetMaker/StJetSimuUtil/StJetSimuReader.h"
00029
00030
00031 #include "StSpinPool/StJets/StJet.h"
00032 #include "StSpinPool/StJets/StJets.h"
00033 #include "StSpinPool/StJetSkimEvent/StJetSkimEvent.h"
00034
00035 ClassImp(StJetSimuReader)
00036
00037
00038 bool verifyJet(StJets* stjets, int ijet);
00039
00040 void dumpProtojetToStream(int event, ostream& os, StJets* stjets);
00041
00042 string idString(TrackToJetIndex* t);
00043
00044
00045 StJetSimuReader::StJetSimuReader(const char* name, StMuDstMaker* uDstMaker)
00046 : StMaker(name), mFile(0), mTree(0), mSkimFile(0), mSkimTree(0),
00047 mDstMaker(uDstMaker), mCounter(0), mValid(false), mSkimEvent(0), mOfstream(0)
00048 {
00049 LOG_DEBUG <<"StJetSimuReader::StJetSimuReader()"<<endm;
00050
00051 }
00052
00053 StJetSimuReader::~StJetSimuReader()
00054 {
00055 LOG_DEBUG <<"StJetSimuReader::~StJetSimuReader()"<<endm;
00056 }
00057
00058
00059 Int_t StJetSimuReader::Init()
00060 {
00061 LOG_DEBUG <<"StJetSimuReader::Init()"<<endm;
00062 return kStOk;
00063 }
00064
00065 void StJetSimuReader::InitFile(const char* file)
00066 {
00067 LOG_DEBUG <<"StJetSimuReader::InitFile()"<<endm;
00068
00069 LOG_DEBUG <<"open file:\t"<<file<<"\tfor reading"<<endm;
00070 mFile = new TFile(file,"READ");
00071 assert(mFile);
00072
00073 LOG_DEBUG <<"recover jet tree"<<endm;
00074 TObject* tree = mFile->Get("jet");
00075 TTree* t = dynamic_cast<TTree*>(tree);
00076 assert(t);
00077 mTree = t;
00078
00079 LOG_DEBUG <<"\tset tree pointer"<<endm;
00080 LOG_DEBUG <<"Number of entries in jet tree:\t"<<t->GetEntries();
00081
00082 LOG_DEBUG <<"\tGet Jet Branches"<<endm;
00083 TObjArray* branches = t->GetListOfBranches();
00084 if (!branches) {LOG_DEBUG <<"StJetSimuReader::InitFile(). Null branches"<<endm; abort();}
00085
00086 LOG_DEBUG <<"\tLoop on branches"<<endm;
00087
00088 for (int i=0; i<branches->GetLast()+1; ++i) {
00089 TBranch* branch = dynamic_cast<TBranch*>((*branches)[i]);
00090 if (!branch) {LOG_DEBUG <<"StJetSimuReader::InitFile(). Null branch"<<endm; abort();}
00091 string bname( branch->GetName() );
00092 LOG_DEBUG <<"\t--- Found branch:\t"<<bname<<endm;
00093
00094 if ( (bname.find("Cone")!=bname.npos) || (bname.find("cone")!=bname.npos) ) {
00095 LOG_DEBUG <<"\t\tcreate StJets object for branch:\t"<<bname<<endm;
00096
00097
00098 StJets* jets = new StJets();
00099 jets->Clear();
00100 mStJetsMap[bname] = jets;
00101 LOG_DEBUG <<"\t\tset branch address for branch:\t"<<bname.c_str()<<endm;
00102 t->SetBranchStatus(bname.c_str(), 1);
00103 t->SetBranchAddress(bname.c_str(), &jets);
00104 branch->SetAutoDelete(true);
00105 }
00106 }
00107
00108 if (0) {
00109 string jetCheck(file);
00110 jetCheck += ".read.txt";
00111 mOfstream = new ofstream(jetCheck.c_str());
00112 }
00113
00114 LOG_DEBUG <<"\tfinished!"<<endm;
00115
00116 return ;
00117 }
00118
00119 void StJetSimuReader::InitJetSkimFile(const char* file)
00120 {
00121 LOG_DEBUG <<" StJetSimuReader::InitJetSkimFile(const char* file)"<<endm;
00122 LOG_DEBUG <<"open file:\t"<<file<<"\tfor reading"<<endm;
00123
00124 mSkimFile = new TFile(file,"READ");
00125 assert(mSkimFile);
00126
00127 TObject* t = mSkimFile->Get("jetSkimTree");
00128 assert(t);
00129
00130 mSkimTree = dynamic_cast<TTree*>(t);
00131 assert(mSkimTree);
00132
00133 mSkimEvent = new StJetSkimEvent();
00134
00135
00136 mSkimTree->SetBranchStatus("skimEventBranch", 1);
00137 mSkimTree->SetBranchAddress("skimEventBranch", &mSkimEvent);
00138
00139 LOG_DEBUG <<"successfully recovered tree from file."<<endm;
00140 }
00141
00142
00143 void StJetSimuReader::InitTree(TTree* tree)
00144 {
00145 LOG_DEBUG <<"StJetSimuReader::InitTree()"<<endm;
00146 LOG_DEBUG <<"\tset tree pointer"<<endm;
00147 LOG_DEBUG <<"Number of entries in tree:\t"<<tree->GetEntries();
00148
00149 TList* friendmist = tree->GetListOfFriends();
00150 TTree* t=0;
00151 for (int j=0; j<friendmist->GetSize()+1; ++j) {
00152 TFriendElement* fr = static_cast<TFriendElement*>( friendmist->At(j) );
00153 string tree_name( fr->GetTreeName() );
00154 if (tree_name == "jet") {
00155 t = fr->GetTree();
00156 break;
00157 }
00158 }
00159 assert(t);
00160
00161 LOG_DEBUG <<"\tGet Branches"<<endm;
00162 TObjArray* branches = t->GetListOfBranches();
00163 if (!branches) {LOG_DEBUG <<"StJetSimuReader::InitFile(). Null branches"<<endm; abort();}
00164
00165 LOG_DEBUG <<"\tLoop on branches"<<endm;
00166
00167 for (int i=0; i<branches->GetLast()+1; ++i) {
00168 TBranch* branch = dynamic_cast<TBranch*>((*branches)[i]);
00169 if (!branch) {LOG_DEBUG <<"StJetSimuReader::InitFile(). Null branch"<<endm; abort();}
00170 string bname( branch->GetName() );
00171 LOG_DEBUG <<"\t--- Found branch:\t"<<bname<<endm;
00172
00173 if ( (bname.find("jet")!=bname.npos) || (bname.find("Jet")!=bname.npos) ) {
00174 LOG_DEBUG <<"\t\tcreate StJets object for branch:\t"<<bname<<endm;
00175
00176
00177 StJets* jets = new StJets();
00178 jets->Clear();
00179 mStJetsMap[bname] = jets;
00180 LOG_DEBUG <<"\t\tset branch address for branch:\t"<<bname.c_str()<<endm;
00181 t->SetBranchStatus(bname.c_str(), 1);
00182 t->SetBranchAddress(bname.c_str(), &jets);
00183 }
00184 }
00185
00186 LOG_DEBUG <<"\tfinished!"<<endm;
00187
00188 return ;
00189 }
00190
00191
00192 int StJetSimuReader::preparedForDualRead()
00193 {
00194 LOG_DEBUG <<"StJetSimuReader::preparedForDualRead()"<<endm;
00195 int jetStatus = (mTree) ? 1 : 0;
00196 int skimStatus = (mSkimTree) ? 1 : 0;
00197 LOG_DEBUG <<"jet tree status:\t"<<jetStatus<<endm;
00198 LOG_DEBUG <<"skim tree status:\t"<<skimStatus<<endm;
00199 assert(mTree);
00200 assert(mSkimTree);
00201
00202 int nJetEvents = mTree->GetEntries();
00203 int nSkimEvents = mSkimTree->GetEntries();
00204 LOG_DEBUG <<"nJetEvents:\t"<<nJetEvents<<endm;
00205 LOG_DEBUG <<"nSkimEvents:\t"<<nSkimEvents<<endm;
00206 assert(nJetEvents==nSkimEvents);
00207 assert(mSkimEvent);
00208
00209 LOG_DEBUG <<"Congratulations, you are ready to process events!"<<endm;
00210 mValid = true;
00211 return 1;
00212 }
00213
00214
00215 Int_t StJetSimuReader::Make()
00216 {
00217 int status = mTree->GetEntry(mCounter);
00218 if (status<0) {
00219 LOG_DEBUG <<"StJetSimuReader::getEvent(). ERROR:\tstatus<0. return null"<<endm;
00220 }
00221
00222 if (mSkimTree) {
00223 status = mSkimTree->GetEntry(mCounter);
00224 }
00225
00226 mCounter++;
00227
00228 return kStOk;
00229 }
00230
00231 Int_t StJetSimuReader::Finish()
00232 {
00233 if (mOfstream) {
00234 delete mOfstream;
00235 mOfstream=0;
00236 }
00237
00238 return kStOk;
00239 }
00240
00241 void StJetSimuReader::exampleSimuAna()
00242 {
00243 LOG_DEBUG <<"StJetSimuReader::exampleSimuAna()"<<endm;
00244
00245 StJetSkimEvent* skEv = skimEvent();
00246
00247
00248 LOG_DEBUG <<"fill/run/event:\t"<<skEv->fill()<<"\t"<<skEv->runId()<<"\t"<<skEv->eventId()<<endm;
00249 LOG_DEBUG <<"fileName:\t"<<skEv->mudstFileName().GetString()<<endm;
00250
00251
00252 LOG_DEBUG <<"Ebbc:\t"<<skEv->eBbc()<<"\tWbbc:\t"<<skEv->wBbc()<<"\tbbcTimeBin:\t"<<skEv->bbcTimeBin()<<endm;
00253
00254
00255 const float* bestPos = skEv->bestVert()->position();
00256 LOG_DEBUG <<"best Vertex (x,y,z):\t"<<bestPos[0]<<"\t"<<bestPos[1]<<"\t"<<bestPos[2]<<endm;
00257
00258
00259 const StPythiaEvent *pyEv=skEv->mcEvent();
00260 LOG_DEBUG <<" subprocess Id="<<pyEv->processId()<<endm;
00261 LOG_DEBUG <<" s =" <<pyEv->s()<<endm;
00262 LOG_DEBUG <<" t =" <<pyEv->t()<<endm;
00263 LOG_DEBUG <<" u =" <<pyEv->u()<<endm;
00264 LOG_DEBUG <<" x1 =" <<pyEv->x1()<<endm;
00265 LOG_DEBUG <<" x2 =" <<pyEv->x2()<<endm;
00266 LOG_DEBUG <<" Q2 =" <<pyEv->Q2()<<endm;
00267 LOG_DEBUG <<" Partonic pT="<<pyEv->pt()<<endm;
00268 LOG_DEBUG <<" Partonic cos(th)="<<pyEv->cosTheta()<<endm;
00269 LOG_DEBUG <<" Partonic Asym="<<pyEv->partonALL()<<endm;
00270 LOG_DEBUG <<" Polarzied PDF for x1="<<pyEv->dF1(StPythiaEvent::STD)<<endm;
00271 LOG_DEBUG <<" Polarzied PDF for x2="<<pyEv->dF2(StPythiaEvent::STD)<<endm;
00272 LOG_DEBUG <<" UnPolarzied PDF for x1="<<pyEv->f1(StPythiaEvent::STD)<<endm;
00273 LOG_DEBUG <<" UnPolarzied PDF for x2="<<pyEv->f1(StPythiaEvent::STD)<<endm;
00274 LOG_DEBUG <<" Weighting =(dF1*dF2*partonALL)/(f1*f2)="<<pyEv->ALL(StPythiaEvent::STD)<<endm;
00275
00276
00277 const TClonesArray* trigs = skEv->triggers();
00278 assert(trigs);
00279 int nTrigs = trigs->GetLast()+1;
00280 for (int i=0; i<nTrigs; ++i) {
00281 StJetSkimTrig* trig = static_cast<StJetSkimTrig*>( (*trigs)[i] );
00282 assert(trig);
00283 if (trig->shouldFire() != 0) {
00284 LOG_DEBUG <<"\tTriggerd:\t"<<trig->trigId()<<endm;
00285 }
00286 }
00287
00288
00289 for (JetBranchesMap::iterator it=mStJetsMap.begin(); it!=mStJetsMap.end(); ++it) {
00290
00291 LOG_DEBUG <<"Found\t"<<(*it).first<<endm;
00292 if ((*it).first!="PythiaConeR04") continue;
00293
00294 StJets* stjets = (*it).second;
00295 int nJets = stjets->nJets();
00296
00297
00298 assert(stjets->runId()==skEv->runId());
00299 assert(stjets->eventId()==skEv->eventId());
00300
00301 TClonesArray* jets = stjets->jets();
00302
00303 for (int ijet=0;ijet<nJets;ijet++){
00304 StJet* j = static_cast<StJet*>( (*jets)[ijet] );
00305 assert(j);
00306 assert(verifyJet(stjets, ijet));
00307 LOG_DEBUG<<"%PYpT="<<j->Pt()<<" PYeta="<<j->Eta()<<" PYphi="<<j->Phi()<<endm;
00308 }
00309
00310 }
00311
00312
00313 for (JetBranchesMap::iterator it=mStJetsMap.begin(); it!=mStJetsMap.end(); ++it) {
00314
00315 LOG_DEBUG <<"Found\t"<<(*it).first<<endm;
00316 if ((*it).first!="MkConeR04") continue;
00317
00318
00319 StJets* stjets = (*it).second;
00320 int nJets = stjets->nJets();
00321 LOG_DEBUG <<"Found\t"<<nJets<<"\tjets from:\t"<<(*it).first<<endm;
00322
00323
00324 assert(stjets->runId()==skEv->runId());
00325 assert(stjets->eventId()==skEv->eventId());
00326 TClonesArray* jets = stjets->jets();
00327
00328 for(int ijet=0; ijet<nJets; ++ijet){
00329 StJet* j = static_cast<StJet*>( (*jets)[ijet] );
00330 assert(j);
00331 assert(verifyJet(stjets, ijet));
00332 LOG_DEBUG <<"jet:\t"<<ijet<<"\tEjet:\t"<<j->E()<<"\tEta:\t"<<j->Eta()<<"\tPhi:\t"<<j->Phi()<<endm;
00333
00334
00335 typedef vector<TLorentzVector*> TrackToJetVec;
00336 TrackToJetVec particles = stjets->particles(ijet);
00337 for (TrackToJetVec::iterator it=particles.begin(); it!=particles.end(); ++it) {
00338 TLorentzVector* t2j = (*it);
00339 assert(t2j);
00340 if (dynamic_cast<TrackToJetIndex*>(t2j)) {LOG_DEBUG<<"TPC track pT="<<t2j->Pt()<<endm;}
00341 if (dynamic_cast<TowerToJetIndex*>(t2j)) {LOG_DEBUG<<"TOWER track eT="<<t2j->Et()<<endm;}
00342 }
00343 }
00344 }
00345 }
00346
00347
00348
00349
00350
00351
