00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 #include <iostream>
00036 #include "StEventTypes.h"
00037 #include "Stypes.h"
00038 #include "StThreeVectorF.hh"
00039 #include "StHelix.hh"
00040 #include "StTrackGeometry.h"
00041 #include "StEventUtilities/StuRefMult.hh"
00042 #include "PhysicalConstants.h"
00043 #include "StPhysicalHelixD.hh"
00044 #include "StTofUtil/StTofGeometry.h"
00045 #include "TNtuple.h"
00046 #include "TFile.h"
00047 #include "StMemoryInfo.hh"
00048 #include "StTimer.hh"
00049 #include "StTofUtil/tofPathLength.hh"
00050 #include "StTofpNtupleMaker.h"
00051
00052
00053 ClassImp(StTofpNtupleMaker)
00054
00055
00057 StTofpNtupleMaker::StTofpNtupleMaker(const Char_t *name){
00058 mTupleFileName="tofntuple.root";
00059 setValidAdcRange(30,1200);
00060 setValidTdcRange(1,2047);
00061 doPrintMemoryInfo = kFALSE;
00062 doPrintCpuInfo = kFALSE;
00063 }
00064
00066 StTofpNtupleMaker::~StTofpNtupleMaker(){ }
00067
00068
00069
00071 Int_t StTofpNtupleMaker::Init(){
00072
00073 if (mTupleFileName!="") bookNtuples();
00074
00075 mTofGeom = new StTofGeometry();
00076 mTofGeom->initDaqMap();
00077
00078 mAcceptedEvents = 0;
00079 mPvpdEntries = 0;
00080 mTofpEvents = 0;
00081 mTofpEntries = 0;
00082
00083 return kStOK;
00084 }
00085
00087 Int_t StTofpNtupleMaker::Finish(){
00088 if (!(mTupleFileName=="")){
00089 mTupleFile->Write();
00090 mTupleFile->Close();
00091 LOG_INFO << "StTofpNtupleMaker::Finish() ntuple file "
00092 << mTupleFileName << " closed." << endm;
00093 }
00094
00095 LOG_INFO << "StTofpNtupleMaker -- statistics" << endm;
00096 LOG_INFO << " accepted events : " << mAcceptedEvents << endm;
00097 LOG_INFO << " pVPD entries : " << mPvpdEntries << endm;
00098 LOG_INFO << " TOFp entries/events : " << mTofpEntries << "/" << mTofpEvents << endm;
00099 return kStOK;
00100 }
00101
00102
00103
00105 Int_t StTofpNtupleMaker::Make(){
00106 LOG_INFO << "StTofpNtupleMaker -- welcome" << endm;
00107
00108 StEvent *event = (StEvent *) GetInputDS("StEvent");
00109
00110
00111
00112 if (!event ||
00113 !event->tofCollection() ||
00114 !event->tofCollection()->dataPresent() ||
00115 !event->primaryVertex()){
00116 LOG_INFO << "StTofpNtupleMaker -- nothing to do ... bye-bye"<< endm;
00117 return kStOK;
00118 }
00119
00120 mAcceptedEvents++;
00121
00122 StTimer timer;
00123 if (doPrintCpuInfo) timer.start();
00124 if (doPrintMemoryInfo) StMemoryInfo::instance()->snapshot();
00125
00126
00127 mYear2= (event->runId()<4000000);
00128 mYear3= (event->runId()>4000000&&event->runId()<5000000);
00129 mYear4= (event->runId()>5000000);
00130
00131
00132
00133
00134
00135
00136 float xvtx = event->primaryVertex()->position().x();
00137 float yvtx = event->primaryVertex()->position().y();
00138 float zvtx = event->primaryVertex()->position().z();
00139 StL0Trigger* pTrigger = event->l0Trigger();
00140 unsigned int triggerWord = 0;
00141 if (pTrigger) triggerWord = pTrigger->triggerWord();
00142
00143
00144 StTriggerDetectorCollection *theTriggers = event->triggerDetectorCollection();
00145 float zdcSumEast(-1.), zdcSumWest(-1), ctbSum(-1);
00146 if (theTriggers){
00147 StCtbTriggerDetector &theCtb = theTriggers->ctb();
00148 StZdcTriggerDetector &theZdc = theTriggers->zdc();
00149 zdcSumEast = theZdc.adcSum(east);
00150 zdcSumWest = theZdc.adcSum(west);
00151 ctbSum = 0;
00152 for (unsigned int islat=0; islat<theCtb.numberOfSlats(); islat++)
00153 for (unsigned int itray=0; itray<theCtb.numberOfTrays(); itray++)
00154 ctbSum += theCtb.mips(itray, islat, 0);
00155 }
00156
00157
00158
00159 int refmult(0);
00160 bool richTofMuDST = (event->summary()->numberOfExoticTracks() == -999);
00161 if (richTofMuDST)
00162 refmult = event->summary()->numberOfGoodTracks();
00163 else
00164 refmult = uncorrectedNumberOfPrimaries(*event);
00165 if (Debug()){
00166 LOG_INFO << " #Tracks :" << event->summary()->numberOfTracks()
00167 << "\n #goodPrimaryTracks:" << event->summary()->numberOfGoodPrimaryTracks()
00168 << "\n #uncorr.prim.tracks :" << refmult << endm;
00169 if (!richTofMuDST)
00170 LOG_INFO << " #goodTracks (global):" << event->summary()->numberOfGoodTracks() << endm;
00171 }
00172
00173
00174 StTofCollection *theTof = event->tofCollection();
00175 getTofData(theTof);
00176
00177
00178 float sumAdcTofp(0.); int nAdcTofp(0), nTdcTofp(0), nAdcTdcTofp(0);
00179 for (int i=0;i<NTOFP;i++){
00180 sumAdcTofp+=mTofpAdc[i];
00181 bool tdcValid=validTdc(mTofpTdc[i]);
00182 bool adcValid=validAdc(mTofpAdc[i]);
00183 if (tdcValid) nTdcTofp++;
00184 if (adcValid) nAdcTofp++;
00185 if (adcValid && tdcValid) nAdcTdcTofp++;
00186 }
00187 if (Debug())
00188 LOG_INFO << " TOFp #Adc:" << nAdcTofp << " #Tdc:" << nTdcTofp << endm;
00189
00190
00191
00192
00193
00194 if (!(mTupleFileName=="")){
00195 int k(0);
00196 float tuple[31];
00197 tuple[k++] = event->runId();
00198 tuple[k++] = event->id();
00199 tuple[k++] = triggerWord;
00200 tuple[k++] = event->summary()->magneticField();
00201 tuple[k++] = zvtx;
00202 tuple[k++] = event->primaryVertex()->chiSquared();
00203 tuple[k++] = ctbSum;
00204 tuple[k++] = zdcSumEast;
00205 tuple[k++] = zdcSumWest;
00206 tuple[k++] = refmult;
00207 tuple[k++] = event->summary()->numberOfGoodPrimaryTracks();
00208 tuple[k++] = sumAdcTofp;
00209 tuple[k++] = nTdcTofp;
00210 for (int i=0;i<NPVPD;i++) tuple[k++] = mPvpdTdc[i];
00211 for (int i=0;i<NPVPD;i++) tuple[k++] = mPvpdAdc[i];
00212 for (int i=0;i<NPVPD;i++) tuple[k++] = mPvpdAdcLoRes[i];
00213
00214 LOG_INFO << " pVPD update ..." << endm;
00215 mPvpdTuple->Fill(tuple);
00216 mPvpdEntries++;
00217 }
00218
00219
00220
00221
00222
00223
00224
00225 if (event->tofCollection()->slatsPresent()){
00226
00227 mTofpEvents++;
00228 int entriesThisEvent(0);
00229
00230
00231 StSPtrVecTofSlat& slatTofVec = theTof->tofSlats();
00232 for (size_t i = 0; i < slatTofVec.size(); i++) {
00233 StTofSlat *thisSlat = slatTofVec[i];
00234 StTrack *thisTrack = thisSlat->associatedTrack();
00235 StTrackGeometry *theTrackGeometry =
00236 (mOuterTrackGeometry)?thisTrack->outerGeometry():thisTrack->geometry();
00237
00238
00239 double pathLength = tofPathLength(&event->primaryVertex()->position(),
00240 &thisSlat->position(),
00241 theTrackGeometry->helix().curvature());
00242 const StThreeVectorF momentum = theTrackGeometry->momentum();
00243
00244
00245 float dedx(0.), cherang(0);
00246 int dedx_np(0), cherang_nph(0);
00247 StSPtrVecTrackPidTraits& traits = thisTrack->pidTraits();
00248 for (unsigned int it=0;it<traits.size();it++){
00249 if (traits[it]->detector() == kTpcId){
00250 StDedxPidTraits *pid = dynamic_cast<StDedxPidTraits*>(traits[it]);
00251 if (pid && pid->method() ==kTruncatedMeanId){
00252 dedx = pid->mean();
00253 dedx_np = pid->numberOfPoints();
00254 }
00255 } else if (traits[it]->detector() == kRichId){
00256 StRichPidTraits *pid = dynamic_cast<StRichPidTraits*>(traits[it]);
00257 if (pid){
00258 StRichSpectra* pidinfo = pid->getRichSpectra();
00259 if (pidinfo && pidinfo->getCherenkovPhotons()>2){
00260 cherang = pidinfo->getCherenkovAngle();
00261 cherang_nph = pidinfo->getCherenkovPhotons();
00262 }
00263 }
00264 }
00265 }
00266
00267
00268 mSlatData.run = event->runId();
00269 mSlatData.evt = event->id();
00270 mSlatData.trgword = triggerWord;
00271 mSlatData.magfield = event->summary()->magneticField();
00272 mSlatData.ctbsum = ctbSum;
00273 mSlatData.zdcsum = zdcSumEast + zdcSumWest;
00274 mSlatData.xvtx = xvtx;
00275 mSlatData.yvtx = yvtx;
00276 mSlatData.zvtx = zvtx;
00277 mSlatData.zvtxchi2 = event->primaryVertex()->chiSquared();
00278 mSlatData.refmult = refmult;
00279 mSlatData.nprimary = event->summary()->numberOfGoodPrimaryTracks();
00280 mSlatData.meanpt = event->summary()->meanPt();
00281 mSlatData.te1 = (int)mPvpdTdc[0]; mSlatData.te2 = (int)mPvpdTdc[1];
00282 mSlatData.te3 = (int)mPvpdTdc[2]; mSlatData.tw1 = (int)mPvpdTdc[3];
00283 mSlatData.tw2 = (int)mPvpdTdc[4]; mSlatData.tw3 = (int)mPvpdTdc[5];
00284 if (mYear3||mYear4) {
00285 mSlatData.ae1 = (int)mPvpdAdcLoRes[0]; mSlatData.ae2 = (int)mPvpdAdcLoRes[1];
00286 mSlatData.ae3 = (int)mPvpdAdcLoRes[2]; mSlatData.aw1 = (int)mPvpdAdcLoRes[3];
00287 mSlatData.aw2 = (int)mPvpdAdcLoRes[4]; mSlatData.aw3 = (int)mPvpdAdcLoRes[5];
00288 } else {
00289 mSlatData.ae1 = (int)mPvpdAdc[0]; mSlatData.ae2 = (int)mPvpdAdc[1];
00290 mSlatData.ae3 = (int)mPvpdAdc[2]; mSlatData.aw1 = (int)mPvpdAdc[3];
00291 mSlatData.aw2 = (int)mPvpdAdc[4]; mSlatData.aw3 = (int)mPvpdAdc[5];
00292 }
00293 mSlatData.slat = thisSlat->slatIndex();
00294 mSlatData.tdc = thisSlat->tdc();
00295 mSlatData.adc = thisSlat->adc();
00296 mSlatData.hitprof = thisSlat->hitProf();
00297 mSlatData.matchflag = thisSlat->matchFlag();
00298 mSlatData.zhit = thisSlat->zHit();
00299 mSlatData.trackId = (Int_t)thisTrack->key();
00300 mSlatData.ntrackpoints= thisTrack->detectorInfo()->numberOfPoints(kTpcId);
00301 mSlatData.nfitpoints = thisTrack->fitTraits().numberOfFitPoints(kTpcId);
00302 mSlatData.r_last = thisTrack->detectorInfo()->lastPoint().perp();
00303 mSlatData.chi2 = thisTrack->fitTraits().chi2(0);
00304
00305 mSlatData.s = fabs(pathLength);
00306 mSlatData.p = momentum.mag()* theTrackGeometry->charge();
00307 mSlatData.pt = momentum.perp();
00308 mSlatData.px = momentum.x();
00309 mSlatData.py = momentum.y();
00310 mSlatData.pz = momentum.z();
00311 mSlatData.eta = momentum.pseudoRapidity();
00312 mSlatData.dedx = dedx;
00313 mSlatData.dedx_np = dedx_np;
00314 mSlatData.cherang = cherang;
00315 mSlatData.cherang_nph = cherang_nph;
00316
00317 mSlatTuple->Fill();
00318 mTofpEntries++;
00319 entriesThisEvent++;
00320 }
00321
00322 LOG_INFO << " TOFp update: " << entriesThisEvent << " entries" <<endm;
00323 }
00324
00325
00326
00327 if (doPrintMemoryInfo) {
00328 StMemoryInfo::instance()->snapshot();
00329 StMemoryInfo::instance()->print();
00330 }
00331 if (doPrintCpuInfo) {
00332 timer.stop();
00333 LOG_INFO << "CPU time for StEventMaker::Make(): "
00334 << timer.elapsedTime() << " sec\n" << endm;
00335 }
00336
00337 LOG_INFO << "StTofpNtupleMaker -- bye-bye" << endm;
00338 return kStOK;
00339 }
00340
00341
00342
00343
00345 Int_t StTofpNtupleMaker::getTofData(StTofCollection* tofCollection){
00346 if (!tofCollection) return kStERR;
00347 StSPtrVecTofData &tofData = tofCollection->tofData();
00348
00349
00350 bool dataOK(true);
00351 if (Debug()) LOG_INFO << "TOF raw data consistency test ..."<< endm;
00352 for (int i=0;i<48;i++){
00353 if (tofData[i]->dataIndex() != mTofGeom->daqToSlatId(i)) {
00354 dataOK = false;
00355 LOG_INFO << "getTofData===>WARNING: " << tofData[i]->dataIndex() << " " << mTofGeom->daqToSlatId(i) << endm;
00356 }
00357
00358 }
00359
00360 for (int i=0;i<NTOFP;i++){
00361 mTofpAdc[i] = tofData[i]->adc();
00362 mTofpTdc[i] = tofData[i]->tdc();
00363 }
00364
00365 if (mYear2){
00366
00367 float tmp = mTofpAdc[3];
00368 mTofpAdc[3] = mTofpAdc[2];
00369 mTofpAdc[2] = tmp;
00370 }
00371
00372 for (int i=0;i<NPVPD;i++){
00373 mPvpdAdc[i] = tofData[42+i]->adc();
00374 mPvpdTdc[i] = tofData[42+i]->tdc();
00375 if (mYear3||mYear4)
00376 mPvpdAdcLoRes[i] = tofData[54+i]->adc();
00377 }
00378
00379 if (!dataOK) return kStWarn;
00380
00381 return kStOK;
00382 }
00383
00384
00386 void StTofpNtupleMaker::bookNtuples(){
00387 mTupleFile = new TFile(mTupleFileName.c_str(), "RECREATE");
00388 LOG_INFO << "StTofpNtupleMaker::bookNtuples() file "
00389 << mTupleFileName << " opened" << endm;
00390
00391
00392 string varList = "run:evt:trgwrd:magfield:zvtx:zvtxchi2:ctbsum"
00393 ":zdceast:zdcwest:refmult:nprimary:tofpsum"
00394 ":ntofp:te1:te2:te3:tw1:tw2:tw3:ael1:ael2:ael3"
00395 ":awl1:awl1:awl3:ae1:ae2:ae3:aw1:aw2:aw3";
00396 mPvpdTuple = new TNtuple("pvpd","tofp timing",varList.c_str());
00397
00398
00399 mSlatTuple = new TTree("tofp","TOFp slat data");
00400 mSlatTuple->Branch("run",&mSlatData.run,"run/I");
00401 mSlatTuple->Branch("evt",&mSlatData.evt,"evt/I");
00402 mSlatTuple->Branch("trgword",&mSlatData.trgword,"trgword/I");
00403 mSlatTuple->Branch("magfield",&mSlatData.magfield,"magfield/F");
00404 mSlatTuple->Branch("ctbsum",&mSlatData.ctbsum,"ctbsum/F");
00405 mSlatTuple->Branch("zdcsum",&mSlatData.zdcsum,"zdcsum/F");
00406 mSlatTuple->Branch("xvtx",&mSlatData.xvtx,"xvtx/F");
00407 mSlatTuple->Branch("yvtx",&mSlatData.yvtx,"yvtx/F");
00408 mSlatTuple->Branch("zvtx",&mSlatData.zvtx,"zvtx/F");
00409 mSlatTuple->Branch("zvtxchi2",&mSlatData.zvtxchi2,"zvtx/F");
00410 mSlatTuple->Branch("refmult",&mSlatData.refmult,"refmult/I");
00411 mSlatTuple->Branch("nprimary",&mSlatData.nprimary,"nprimary/I");
00412 mSlatTuple->Branch("meanpt",&mSlatData.meanpt,"meanpt/F");
00413 mSlatTuple->Branch("tdcstart",&mSlatData.tdcstart,"tdcstart/F");
00414 mSlatTuple->Branch("pvpd",&mSlatData.te1,"te1/I:te2:te3:tw1:tw2:tw3:ae1:ae2:ae3:aw1:aw2:aw3");
00415 mSlatTuple->Branch("slat",&mSlatData.slat,"slat/I:tdc:adc:hitprof:matchflag:zhit/F:zhitinner/F:zhitouter/F:ss:theta_xy:theta_zr");
00416
00417 mSlatTuple->Branch("track",&mSlatData.trackId,"trackId/I:ntrackpoints/I:nfitpoints:r_last/F:chi2:s:p:pt:px:py:pz:eta:dedx/F:dedx_np/I:cherangle/F:cherangle_nph/I");
00418
00419
00420
00421 mMatchTuple = new TTree("match","TOFp match data");
00422 mMatchTuple->Branch("id",&mMatchData.daqid,"daqid/I");
00423 mMatchTuple->Branch("slat",&mMatchData.nneighbors,"nneighbors/I:zlocal/F:philocal/F:hitprof/I");
00424 mMatchTuple->Branch("track",&mMatchData.nfitpoints,"nfitpoints/I:ntrackpoints/I:maxpoints/I:r_last/F:p/F:pt/F");
00425
00426 mNoMatchTuple = new TTree("nomatch","TOFp no-match data");
00427 mNoMatchTuple->Branch("id",&mMatchData.daqid,"daqid/I");
00428 mNoMatchTuple->Branch("slat",&mMatchData.nneighbors,"nneighbors/I:zlocal/F:philocal/F:hitprof/I");
00429 mNoMatchTuple->Branch("track",&mMatchData.nfitpoints,"nfitpoints/I:ntrackpoints/I:maxpoints/I:r_last/F:p/F:pt/F");
00430
00431 return;
00432 }
00433
00434
00435
