00001
00002
00003
00004
00005
00006
00007 #include "StChain.h"
00008
00009 #include "StEvent/StEventTypes.h"
00010
00011 #include "StEventUtilities/StuRefMult.hh"
00012 #include "StEventUtilities/StuProbabilityPidAlgorithm.h"
00013
00014 #include "StarClassLibrary/StPhysicalHelixD.hh"
00015 #include "StarClassLibrary/StTimer.hh"
00016
00017 #include "StIOMaker/StIOMaker.h"
00018
00019 #include "StStrangeMuDstMaker/StStrangeMuDstMaker.h"
00020 #include "StStrangeMuDstMaker/StStrangeEvMuDst.hh"
00021 #include "StStrangeMuDstMaker/StV0MuDst.hh"
00022 #include "StStrangeMuDstMaker/StV0Mc.hh"
00023 #include "StStrangeMuDstMaker/StXiMuDst.hh"
00024 #include "StStrangeMuDstMaker/StXiMc.hh"
00025 #include "StStrangeMuDstMaker/StKinkMuDst.hh"
00026 #include "StStrangeMuDstMaker/StKinkMc.hh"
00027
00028 #include "StMuDSTMaker/COMMON/StMuException.hh"
00029 #include "StMuDSTMaker/COMMON/StMuEvent.h"
00030 #include "StMuDSTMaker/COMMON/StMuTrack.h"
00031 #include "StMuDSTMaker/COMMON/StMuDebug.h"
00032 #include "StMuDSTMaker/COMMON/StMuCut.h"
00033 #include "StMuDSTMaker/COMMON/StMuDst.h"
00034
00035 #include "StHbtMuDstReader.h"
00036 #include "Infrastructure/StHbtEvent.hh"
00037 #include "StHbtMaker/Base/StHbtEventCut.h"
00038
00039 #include "TFile.h"
00040 #include "TTree.h"
00041 #include "TClass.h"
00042 #include "TChain.h"
00043 #include "TStreamerInfo.h"
00044 #include "TClonesArray.h"
00045
00046 ClassImp(StHbtMuDstReader)
00047
00048 #if !(ST_NO_NAMESPACES)
00049 using namespace units;
00050 #endif
00051
00052
00053
00054
00055
00056 StHbtMuDstReader::StHbtMuDstReader(int mode, int nameMode, const char* dirName, const char* fileName, const char* filter, int maxFiles) :
00057 mStEvent(0), mStStrangeMuDstMaker(0), mIOMaker(0),
00058 mIoMode((ioMode)mode), mIoNameMode((ioNameMode)nameMode),
00059 mDirName(dirName), mFileName(fileName), mFilter(filter), mMaxFiles(maxFiles),
00060 mTrackType(primary), mReadTracks(1),
00061 mReadV0s(1), mReadXis(1), mReadKinks(1), mFinish(0),
00062 mSplit(99), mCompress(9), mBufferSize(65536*4), mHbtEvent(0)
00063 {
00064 mReaderStatus = 0;
00065 mEventCounter=0;
00066 mStMuDst = new StMuDst();
00068 for ( int i=0; i<__NARRAYS__; i++) {
00069 arrays[i] = 0;
00070 mArrays[i]= clonesArray(arrays[i],StMuArrays::arrayTypes[i],StMuArrays::arraySizes[i],StMuArrays::arrayCounters[i]);
00071 }
00073 for ( int i=0; i<__NSTRANGEARRAYS__; i++) {
00074 strangeArrays[i] = 0;
00075 mStrangeArrays[i]= clonesArray(strangeArrays[i],StMuArrays::strangeArrayTypes[i],StMuArrays::strangeArraySizes[i],StMuArrays::strangeArrayCounters[i]);
00076 }
00077
00078 mStMuDst->set(arrays,strangeArrays);
00079
00080 if (mIoMode==ioRead) openRead();
00081 if (mIoMode==ioWrite) mProbabilityPidAlgorithm = new StuProbabilityPidAlgorithm();
00082
00083 }
00084
00085
00086
00087
00088 StHbtMuDstReader::~StHbtMuDstReader(){
00089 }
00090
00091
00092
00093
00094 void StHbtMuDstReader::clear(){
00095 DEBUGMESSAGE1("");
00096
00097
00099 for ( int i=0; i<__NARRAYS__; i++) {
00100 clear(mArrays[i],StMuArrays::arrayCounters[i]);
00101 }
00102 for ( int i=0; i<__NSTRANGEARRAYS__; i++) {
00103 clear(mStrangeArrays[i],StMuArrays::strangeArrayCounters[i]);
00104 }
00105
00106
00107 }
00108
00109
00110
00111 void StHbtMuDstReader::clear(TClonesArray* t, int& counter){
00112 if (t) t->Clear(""); counter=0;
00113 }
00114
00115
00116
00117 TClonesArray* StHbtMuDstReader::clonesArray(TClonesArray* p, const char* type, int size, int& counter) {
00118 DEBUGMESSAGE1("");
00119 if (!p) {
00120 DEBUGVALUE2(type);
00121 p = new TClonesArray(type, size);
00122 counter=0;
00123 }
00124 if (!p) throw StMuExceptionNullPointer("could not create TClonesArray",__PRETTYF__);
00125 return p;
00126 }
00127
00128
00129
00130 int StHbtMuDstReader::Init(){
00131 DEBUGMESSAGE1("");
00132
00133
00134 return 0;
00135 }
00136
00137
00138
00139
00140 void StHbtMuDstReader::Clear(){
00141 DEBUGMESSAGE1("");
00142 clear();
00143 }
00144
00145
00146
00147 StHbtEvent* StHbtMuDstReader::ReturnHbtEvent(){
00148 DEBUGMESSAGE1("");
00149 StTimer timer;
00150 timer.start();
00151 clear();
00152 try {
00153 if (mIoMode == ioRead) read();
00154 }
00155 catch(StMuExceptionEOF e) {
00156 e.print();
00157 mReaderStatus = 2;
00158 return 0;
00159 }
00160 catch(StMuException e) {
00161 e.print();
00162 mReaderStatus = 1;
00163 return 0;
00164 }
00165 DEBUGVALUE1(timer.elapsedTime());
00166 if (mStMuDst) {
00167 mStMuDst->fixTrackIndices();
00168 mHbtEvent = new StHbtEvent(mStMuDst, mTrackType);
00169 }
00170 if (mEventCut){
00171 if (!(mEventCut->Pass(mHbtEvent))){
00172 delete mHbtEvent;
00173 mHbtEvent = 0;
00174 return 0;
00175 }
00176 }
00177 return mHbtEvent;
00178 }
00179
00180
00181
00182 void StHbtMuDstReader::fill(){
00183 DEBUGMESSAGE1("");
00184
00185 if (!mStEvent) {
00186 DEBUGMESSAGE1("no StEvent");
00187 return;
00188 }
00189
00192 if (mProbabilityPidAlgorithm) delete mProbabilityPidAlgorithm;
00193 mProbabilityPidAlgorithm = new StuProbabilityPidAlgorithm(*mStEvent);
00194 StMuTrack::setProbabilityPidAlgorithm(mProbabilityPidAlgorithm);
00195 StMuTrack::setProbabilityPidCentrality(uncorrectedNumberOfNegativePrimaries(*mStEvent));
00196
00197
00198 try {
00199 fillTrees(mStEvent);
00200 }
00201 catch(StMuException e) {
00202 e.print();
00203 throw e;
00204 }
00205 }
00206
00207
00208
00209 void StHbtMuDstReader::write(){
00210 DEBUGMESSAGE1("");
00211
00212
00213 try {
00214 fill();
00215 }
00216 catch (StMuException e) {
00217 return;
00218 }
00219
00220
00221 string ioMakerFileName;
00222 if (mIOMaker) {
00223 ioMakerFileName = string(mIOMaker->GetFile());
00224 }
00225 else {
00226 ioMakerFileName = mDirName+mFileName;
00227 }
00228 DEBUGVALUE1(ioMakerFileName.c_str());
00229
00230
00231 string theFileName = buildFileName(mDirName,basename(ioMakerFileName),".MuDst.root");
00232 if (theFileName != mCurrentFileName) {
00233 closeWrite();
00234 openWrite(theFileName);
00235 mCurrentFileName = theFileName;
00236 }
00237
00238 DEBUGMESSAGE2("now fill tree");
00239 mTTree->Fill();
00240 DEBUGMESSAGE2("tree filled");
00241
00242 return;
00243 }
00244
00245
00246
00247 void StHbtMuDstReader::Finish() {
00248 if (mFinish) {
00249 for ( int i=0; i<10; i++) {
00250 cout << "why are you calling the Finish() again ???????" << endl;
00251 cout << "are you the stupid chain destructor ???????????" << endl;
00252 }
00253 }
00254 else {
00255 if (mIoMode== ioWrite ) closeWrite();
00256 if (mIoMode== ioRead ) closeRead();
00257 mFinish = true;
00258 }
00259 return;
00260 }
00261
00262
00263
00264 void StHbtMuDstReader::openRead() {
00265 DEBUGVALUE1(mDirName.c_str());
00266 DEBUGVALUE1(mFileName.c_str());
00267 DEBUGVALUE1(mFilter.c_str());
00268
00269 makeChain(mDirName.c_str(), mFilter.c_str(),mMaxFiles);
00270
00271
00272 for ( int i=0; i<__NARRAYS__; i++) {
00273 mChain->SetBranchAddress(StMuArrays::arrayNames[i],&mArrays[i]);
00274 }
00275
00276
00277 for ( int i=0; i<__NSTRANGEARRAYS__; i++) {
00278 mChain->SetBranchAddress(StMuArrays::strangeArrayNames[i],&mStrangeArrays[i]);
00279 }
00280
00281 (void*)mChain->GetBranch("");
00282 mTTree = mChain->GetTree();
00283
00284 mStMuDst->set(mArrays,mStrangeArrays);
00285 }
00286
00287
00288
00289 void StHbtMuDstReader::read(){
00290 DEBUGMESSAGE1("");
00291 if ( !(mEventCounter<mChain->GetEntries()) ) throw StMuExceptionEOF("end of input",__PRETTYF__);
00292 mChain->GetEntry(mEventCounter);
00293 DEBUGVALUE2(mChain->GetCurrentFile()->GetName());
00294 mEventCounter++;
00295
00296 return;
00297 }
00298
00299
00300
00301 void StHbtMuDstReader::closeRead(){
00302 }
00303
00304
00305
00306 void StHbtMuDstReader::openWrite(string fileName) {
00307 DEBUGVALUE1(fileName.c_str());
00308
00309 DEBUGMESSAGE2("now create file");
00310 mCurrentFile = new TFile(fileName.c_str(),"RECREATE","StMuDst");
00311
00312 if (!mCurrentFile) throw StMuExceptionNullPointer("no file openend",__PRETTYF__);
00313
00314 mCurrentFile->SetCompressionLevel(mCompress);
00315
00316
00317 DEBUGMESSAGE2("now create trees and branches");
00318
00319 TBranch* branch;
00320 int bufsize = mBufferSize;
00321 if (mSplit) bufsize /= 4;
00322
00323
00324 mTTree = new TTree("MuDst", "StMuDst",mSplit);
00325 if (!mTTree) throw StMuExceptionNullPointer("can not create tree",__PRETTYF__);
00326 mTTree->SetAutoSave(1000000);
00327 DEBUGMESSAGE("arrays");
00328 for ( int i=0; i<__NARRAYS__; i++) {
00329 DEBUGVALUE2(i);
00330 branch = mTTree->Branch(StMuArrays::arrayNames[i],&mArrays[i], bufsize, mSplit);
00331 }
00332
00333
00334 DEBUGMESSAGE("strange arrays");
00335 for ( int i=0; i<__NSTRANGEARRAYS__; i++) {
00336 DEBUGVALUE2(i);
00337 branch = mTTree->Branch(StMuArrays::strangeArrayNames[i],&mStrangeArrays[i], bufsize, mSplit);
00338 }
00339
00340 mCurrentFileName = fileName;
00341 }
00342
00343
00344
00345 void StHbtMuDstReader::closeWrite(){
00346 if (mTTree) mTTree->AutoSave();
00347 if (mCurrentFile) mCurrentFile->Close();
00348 mTTree = 0;
00349 mCurrentFile = 0;
00350 }
00351
00352
00353
00354 void StHbtMuDstReader::fillTrees(StEvent* ev, StMuCut* cut){
00355 DEBUGMESSAGE1("");
00356
00357 try {
00358 fillEvent(ev);
00359 fillL3AlgorithmInfo(ev);
00360 fillDetectorStates(ev);
00361 }
00362 catch(StMuException e) {
00363 e.print();
00364 }
00365
00366 try {
00367 fillTracks(ev,mTrackFilter);
00368 }
00369 catch(StMuException e) {
00370 e.print();
00371 }
00372
00373 try {
00374 fillL3Tracks(ev, mL3TrackFilter);
00375 }
00376 catch(StMuException e) {
00377 e.print();
00378 }
00379
00380 try {
00381 fillStrange(mStStrangeMuDstMaker);
00382 }
00383 catch(StMuException e) {
00384 e.print();
00385 }
00386 }
00387
00388
00389
00390
00391
00392
00393 void StHbtMuDstReader::fillEvent(StEvent* ev, StMuCut* cut) {
00394 DEBUGMESSAGE1("");
00395 StMuEvent typeOfEvent;
00396 if (!ev) throw StMuExceptionNullPointer("no StEvent",__PRETTYF__);
00397 StTimer timer;
00398 timer.start();
00399 if (!cut || cut->pass(ev)) {
00400 DEBUGMESSAGE3("");
00401 addType(mArrays[muEvent],ev,typeOfEvent);
00402 }
00403 timer.stop();
00404 DEBUGVALUE2(timer.elapsedTime());
00405 }
00406
00407
00408
00409 void StHbtMuDstReader::fillL3AlgorithmInfo(StEvent* ev) {
00410 DEBUGMESSAGE1("");
00411 if ( !ev->l3Trigger() ) return;
00412 if ( !ev->l3Trigger()->l3EventSummary()) return;
00413
00414 StTimer timer;
00415 timer.start();
00416 StL3EventSummary* l3 = ev->l3Trigger()->l3EventSummary();
00417 int n = l3->numberOfAlgorithms();
00418 for (int i=0; i<n; i++) {
00419 if (l3->algorithms()[i]->accept())
00420 addType( mArrays[muAccept], *l3->algorithms()[i] );
00421 else
00422 addType( mArrays[muReject], *l3->algorithms()[i] );
00423 }
00424 timer.stop();
00425 DEBUGVALUE2(timer.elapsedTime());
00426 }
00427
00428
00429
00430 void StHbtMuDstReader::fillTracks(StEvent* ev, StMuCut* cut) {
00431 DEBUGMESSAGE1("");
00432 StTimer timer;
00433 timer.start();
00434
00435 StSPtrVecTrackNode& nodes= ev->trackNodes();
00436 DEBUGVALUE2(nodes.size());
00437 for (StSPtrVecTrackNodeConstIterator iter=nodes.begin(); iter!=nodes.end(); iter++) {
00438 addTrackNode(ev, *iter, cut, mArrays[muGlobal], mArrays[muPrimary], mArrays[muOther], false);
00439 }
00440 timer.stop();
00441 DEBUGVALUE2(timer.elapsedTime());
00442 }
00443
00444
00445
00446 void StHbtMuDstReader::fillL3Tracks(StEvent* ev, StMuCut* cut) {
00447 DEBUGMESSAGE1("");
00448 if (!ev->l3Trigger()) return;
00449
00450 StTimer timer;
00451 timer.start();
00452 StSPtrVecTrackNode& nodes= ev->l3Trigger()->trackNodes();
00453 DEBUGVALUE2(nodes.size());
00454 for (StSPtrVecTrackNodeConstIterator iter=nodes.begin(); iter!=nodes.end(); iter++) {
00455 addTrackNode(ev, *iter, cut, mArrays[muL3], 0, 0, true );
00456 }
00457 timer.stop();
00458 DEBUGVALUE2(timer.elapsedTime());
00459 }
00460
00461
00462
00463 void StHbtMuDstReader::fillDetectorStates(StEvent* ev) {
00464 DEBUGMESSAGE1("");
00465 StTimer timer;
00466 timer.start();
00467 for (int i=0; i<StMuArrays::arraySizes[muState]; i++) {
00468 StDetectorState* state = ev->detectorState((StDetectorId) i);
00469 if (state)
00470 addType( mArrays[muState], ev->detectorState((StDetectorId)i) );
00471 }
00472 timer.stop();
00473 DEBUGVALUE2(timer.elapsedTime());
00474 }
00475
00476
00477
00478 void StHbtMuDstReader::addTrackNode(const StEvent* ev, const StTrackNode* node, StMuCut* cut,
00479 TClonesArray* gTCA, TClonesArray* pTCA, TClonesArray* oTCA, bool l3) {
00480 DEBUGMESSAGE3("");
00481 const StTrack* tr=0;
00482
00484 int index2Global =-1;
00485 if (gTCA) {
00486 tr= node->track(global);
00487 if (tr ) index2Global = addTrack(gTCA, ev, tr, cut, -1, l3);
00488 }
00490 int index;
00491 if (pTCA) {
00492 tr = node->track(primary);
00493 if (tr) index = addTrack(pTCA, ev, tr, cut, index2Global, l3);
00494 }
00496 if (oTCA) {
00497 size_t nEntries = node->entries();
00498 for (size_t j=0; j<nEntries; j++) {
00499 tr = node->track(j);
00500 if (tr && (tr->type()!=global) && (tr->type()!=primary) ) {
00501 index = addTrack(oTCA, ev, tr, cut, index2Global, l3);
00502 }
00503 }
00504 }
00505 }
00506
00507
00508
00509 int StHbtMuDstReader::addTrack(TClonesArray* tca, const StEvent* event, const StTrack* track, StMuCut* cut, int index2Global, bool l3) {
00510 DEBUGMESSAGE3("");
00511 StRichSpectra typeOfStRichSpectra;
00512 int index = -1;
00513 int index2RichSpectra=-1;
00515 int counter = tca->GetEntries();
00516 try{
00517 if (cut && !cut->pass(track)) throw StMuExceptionBadValue("failed track cut",__PRETTYF__);
00518
00519
00520 StRichSpectra* rich = richSpectra(track);
00521 if (rich) {
00522 index2RichSpectra = addType( mArrays[muRich], *rich );
00523 }
00524 new((*tca)[counter]) StMuTrack(event, track, 0, index2Global, index2RichSpectra, l3);
00525 index = counter;
00526 }
00527 catch (StMuException e) {
00528 IFDEBUG3(e.print());
00529 }
00530 return index;
00531 }
00532
00533
00534
00535 StRichSpectra* StHbtMuDstReader::richSpectra(const StTrack* track) {
00536 DEBUGMESSAGE3("");
00537 const StPtrVecTrackPidTraits& traits = track->pidTraits(kRichId);
00538 for (StPtrVecTrackPidTraitsConstIterator traitIter=traits.begin();traitIter!=traits.end();++traitIter) {
00539 StRichPidTraits* pid = dynamic_cast<StRichPidTraits*>(*traitIter);
00540 if (pid) return pid->getRichSpectra();
00541 }
00542 return 0;
00543 }
00544
00545
00546
00547 void StHbtMuDstReader::fillStrange(StStrangeMuDstMaker* maker) {
00548 DEBUGMESSAGE2("");
00550 if (!maker) throw StMuExceptionNullPointer("no StrangeMuDstMaker",__PRETTYF__);
00551
00552 StStrangeEvMuDst ev;
00553 StV0MuDst v0;
00554 StStrangeAssoc assoc;
00555 StXiMuDst xi;
00556 StKinkMuDst kink;
00557 StV0Mc v0Mc;
00558 StXiMc xiMc;
00559 StKinkMc kinkMc;
00560
00561 addType(maker->GetEvClonesArray(), mStrangeArrays[0],ev);
00562 addType(maker->GetEvMcArray(), mStrangeArrays[1],ev);
00563
00564 addType(maker->GetV0ClonesArray(), mStrangeArrays[2],v0);
00565 addType(maker->GetV0McArray(), mStrangeArrays[3],v0Mc);
00566 addType(maker->GetV0AssocArray(), mStrangeArrays[4],assoc);
00567
00568 addType(maker->GetXiClonesArray(), mStrangeArrays[5],xi);
00569 addType(maker->GetXiMcArray(), mStrangeArrays[6],xiMc);
00570 addType(maker->GetXiAssocArray(), mStrangeArrays[7],assoc);
00571
00572 addType(maker->GetKinkClonesArray(),mStrangeArrays[8],kink);
00573 addType(maker->GetKinkMcArray(), mStrangeArrays[9],kinkMc);
00574 addType(maker->GetKinkAssocArray(), mStrangeArrays[10],assoc);
00575
00576 }
00577
00578
00579
00580 template <class T>
00581 void StHbtMuDstReader::addType(TClonesArray* tcaFrom, TClonesArray* tcaTo , T t) {
00582 if (tcaFrom && tcaTo) {
00583 int n = tcaFrom->GetEntries();
00584 int counter = tcaTo->GetEntries();
00585 for (int i=0; i<n;i++) {
00586 new((*tcaTo)[counter++]) T( (T&)*tcaFrom->UncheckedAt(i) );
00587 }
00588 }
00589 }
00590
00591
00592
00593 template <class T>
00594 int StHbtMuDstReader::addType(TClonesArray* tcaTo , T t) {
00595 int counter =-1;
00596 if (tcaTo) {
00597 counter = tcaTo->GetEntries();
00598 new((*tcaTo)[counter]) T( t );
00599 }
00600 return counter;
00601 }
00602
00603
00604
00605 template <class T, class U>
00606 int StHbtMuDstReader::addType(TClonesArray* tcaTo , U u, T t) {
00607 int counter =-1;
00608 if (tcaTo) {
00609 counter = tcaTo->GetEntries();
00610 DEBUGMESSAGE("");
00611 new((*tcaTo)[counter]) T(u);
00612 }
00613 return counter;
00614 }
00615
00616
00617
00618 string StHbtMuDstReader::buildFileName(string dir, string fileName, string extention){
00619 DEBUGMESSAGE1("");
00620 fileName = dir + "/" + fileName + extention;
00621 return fileName;
00622 }
00623
00624
00625
00626 string StHbtMuDstReader::basename(string s){
00627 string name(s);
00628 while ( name.find("/") != string::npos ) {
00629 string::size_type pos = name.find("/");
00630 name.erase(0, pos+1 );
00631 }
00632 string::size_type pos = name.find(".event.root");
00633 if (pos != string::npos) name.erase(pos,pos+11);
00634 return name;
00635 }
00636
00637
00638
00639 void StHbtMuDstReader::makeChain(const char* dir, const char* filter, int maxFiles) {
00640 DEBUGMESSAGE1("");
00641 mChain = new TChain("MuDst");
00642
00643
00644
00645 int fileCount(0);
00646 if(strncmp(dir+strlen(dir)-4,".lis",4)==0){
00647 ifstream tF(dir);
00648 char tFileName[500];
00649 int tNFile =0;
00650 tF >> tFileName;
00651 while ((!tF.eof()) && (tNFile<=maxFiles)){
00652 cout << "Add file " << tFileName << endl;
00653 mChain->Add(tFileName);
00654 tNFile++;
00655 tF >> tFileName;
00656 }
00657 fileCount = tNFile;
00658 }
00659 else{
00660 const char* fileName(0);
00661 void *pDir = gSystem->OpenDirectory(dir);
00662 while((fileName = gSystem->GetDirEntry(pDir))){
00663 if(strcmp(fileName,".")==0 || strcmp(fileName,"..")==0) continue;
00664 if(strcmp(fileName,".event.root")==0 || strcmp(fileName,"..")==0) continue;
00665
00666 if(strstr(fileName,filter) && strstr(fileName,".MuDst.root") ){
00667 char* fullFile = gSystem->ConcatFileName(dir,fileName);
00668
00669 cout << fileCount << " " << fullFile << endl;
00670 mChain->Add(fullFile);
00671
00672 delete fullFile;
00673 if(++fileCount >= maxFiles) break;
00674 }
00675 }
00676 }
00677 DEBUGVALUE2(fileCount);
00678 }
00679
00680 void StHbtMuDstReader::setProbabilityPidFile(const char* file) {
00681 if (mProbabilityPidAlgorithm)
00682 mProbabilityPidAlgorithm->readParametersFromFile(file);
00683 }
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
