00001 #include <TROOT.h>
00002 #include <TTree.h>
00003 #include <TH2F.h>
00004 #include <TVector3.h>
00005 #include <TFile.h>
00006 #include <assert.h>
00007
00008 #include <StEvent.h>
00009 #include <StEventTypes.h>
00010
00011 #include <StMuDSTMaker/COMMON/StMuTrack.h>
00012 #include <StMuDSTMaker/COMMON/StMuDst.h>
00013 #include <StMuDSTMaker/COMMON/StMuEvent.h>
00014 #include <StMuDSTMaker/COMMON/StMuTriggerIdCollection.h>
00015 #include <StMuDSTMaker/COMMON/StMuPrimaryVertex.h>
00016
00017 #include <StEmcUtil/geometry/StEmcGeom.h>
00018 #include <StEmcUtil/projection/StEmcPosition.h>
00019 #include <StEmcUtil/database/StBemcTables.h>
00020 #include <StEmcUtil/others/emcDetectorName.h>
00021 #include <StEmcRawMaker/defines.h>
00022 #include <StEvent/StEnumerations.h>
00023 #include <St_db_Maker/St_db_Maker.h>
00024 #include <StEmcADCtoEMaker/StEmcADCtoEMaker.h>
00025 #include <StDetectorDbMaker/StDetectorDbTriggerID.h>
00026
00027
00028 #include <StEmcPool/StPhotonCommon/MyEvent.h>
00029 #include <StEmcPool/StPhotonCommon/MyPoint.h>
00030 #include <StEmcPool/StPhotonCommon/MyMcTrack.h>
00031
00032 #include "StMyEventMaker.h"
00033
00034 ClassImp(StMyEventMaker)
00035
00036 StMyEventMaker::StMyEventMaker(const char *name,const char *outputFile,
00037 const char *type,const char *coll,Bool_t debug):StMaker(name)
00038 {
00039 mFileName=outputFile;
00040
00041 if(strcmp(type,"real")!=0){
00042 cout<<"not running on real data??"<<endl;
00043 assert(0);
00044 }
00045
00046 mDAU=kFALSE;
00047 mPP04=kFALSE;
00048 mPP05=kFALSE;
00049 mAUAU200=kFALSE;
00050 if(strcmp(coll,"dAu")==0){
00051 mDAU=kTRUE;
00052 mTrig[0]=2001;
00053 mTrig[1]=2003;
00054 mTrig[2]=2201;
00055 mTrig[3]=2202;
00056 }
00057 if(strcmp(coll,"pp04")==0){
00058 mPP04=kTRUE;
00059 mTrig[0]=45010;
00060 mTrig[1]=10;
00061 mTrig[2]=45201;
00062 mTrig[3]=45202;
00063 }
00064 if(strcmp(coll,"pp05")==0){
00065 mPP05=kTRUE;
00066 mTrig[0]=96011;
00067 mTrig[1]=-10000;
00068 mTrig[2]=96201;
00069 mTrig[3]=96211;
00070 }
00071 if(strcmp(coll,"auau200")==0){
00072 mAUAU200=kTRUE;
00073 mTrig[0]=15003;
00074 mTrig[1]=15007;
00075 mTrig[2]=15203;
00076 mTrig[3]=-10000;
00077 }
00078
00079 mDebug=debug;
00080
00081 mBemcTables=new StBemcTables();
00082 mEvent=new MyEvent();
00083 }
00084
00085 StMyEventMaker::~StMyEventMaker()
00086 {
00087
00088 }
00089
00090 Int_t StMyEventMaker::Init()
00091 {
00092
00093 h_EvSum=new TH1F("h_EvSum","see code for details",22,-1.5,20.5);
00094
00095 mEventTree=new TTree("mEventTree","tree with my event structure");
00096 mEventTree->Branch("branch","MyEvent",&mEvent);
00097
00098 mN=0;
00099 mRunId=0;
00100 mRunPrev=0;
00101 mEventId=0;
00102 mTrigger=0;
00103 mDate=0;
00104 mTime=0;
00105 mPs_mb=0;
00106 mPs_mb2=0;
00107 mPs_ht1=0;
00108 mPs_ht2=0;
00109
00110 return StMaker::Init();
00111 }
00112
00113 Int_t StMyEventMaker::Make()
00114 {
00115 mN++;
00116 if(mDebug) cout<<endl<<"#### processing event "<<mN<<endl<<endl;
00117 h_EvSum->Fill(1);
00118
00119 StEmcGeom *bemcGeom=StEmcGeom::getEmcGeom("bemc");
00120
00121 StMuDst* mu=(StMuDst*)GetInputDS("MuDst");
00122 if(!mu){
00123 if(mDebug) cout << "********* no StMuDst*"<<endl;
00124 return kStOK;
00125 }
00126
00127 if(mDebug) cout<<"number of primverts: "<<mu->numberOfPrimaryVertices()<<endl;
00128
00129 StMuEvent *muEvent=(StMuEvent*)mu->event();
00130 if(!muEvent){
00131 if(mDebug) cout << "********* no StMuEvent*"<<endl;
00132 return kStOK;
00133 }
00134
00135 mRunPrev=mRunId;
00136 mRunId=muEvent->runId();
00137 mEventId=muEvent->eventId();
00138 if(mDebug) cout<<"starting run: "<<mRunId<<" and event: "<<mEventId<<endl;
00139
00140 if(mAdcMaker->isCorrupted()){
00141 if(mDebug) cout<<"event corrupted!"<<endl;
00142 h_EvSum->Fill(-1);
00143 return kStOK;
00144 }
00145
00146 if(mRunId!=mRunPrev){
00147 StDetectorDbTriggerID& DetDbTrigId = *(StDetectorDbTriggerID::instance());
00148 for(UInt_t i=0;i<DetDbTrigId.getL0NumRows();i++){
00149 if(mDebug) cout<<"prescales: "<<DetDbTrigId.getPsL0(i)<<endl<<endl;
00150 if(DetDbTrigId.getL0OfflineTrgId(i)==mTrig[0])
00151 mPs_mb=DetDbTrigId.getPsL0(i);
00152 if(DetDbTrigId.getL0OfflineTrgId(i)==mTrig[1])
00153 mPs_mb2=DetDbTrigId.getPsL0(i);
00154 if(DetDbTrigId.getL0OfflineTrgId(i)==mTrig[2])
00155 mPs_ht1=DetDbTrigId.getPsL0(i);
00156 if(DetDbTrigId.getL0OfflineTrgId(i)==mTrig[3])
00157 mPs_ht2=DetDbTrigId.getPsL0(i);
00158 }
00159
00160 mBemcTables->loadTables((StMaker*)this);
00161 }
00162
00163 mDate=(Int_t)mDbMaker->GetDateTime().GetDate();
00164 mTime=(Int_t)mDbMaker->GetDateTime().GetTime();
00165
00166
00167 Float_t mips=-1.;
00168 Float_t zdcEast=0.;
00169 Float_t zdcWest=0.;
00170 Float_t zdcVertex=0.;
00171 Float_t bbcEast=0.;
00172 Float_t bbcWest=0.;
00173 StCtbTriggerDetector &CTB=muEvent->ctbTriggerDetector();
00174 mips=CTB.mips(0);
00175 StZdcTriggerDetector &ZDC=muEvent->zdcTriggerDetector();
00176 zdcEast=ZDC.adcSum(east);
00177 zdcWest=ZDC.adcSum(west);
00178 zdcVertex=ZDC.vertexZ();
00179 StBbcTriggerDetector &BBC=muEvent->bbcTriggerDetector();
00180 bbcEast=BBC.adcSumEastLarge();
00181 bbcWest=BBC.adcSumWestLarge();
00182
00183
00184 int west=BBC.tdcEarliestWest();
00185 int east=BBC.tdcEarliestEast();
00186 int diff=west-east;
00187 float BBCVertexZ = 10.77 + 2.82*diff;
00188
00189
00190
00191
00192
00193
00194 mEvent=new MyEvent(mRunId,mEventId,mDate,mTime);
00195 if(!mEvent){
00196 if(mDebug) cout<<"couldn't create MyEvent*"<<endl;
00197 return kStErr;
00198 }
00199
00200
00201 StThreeVectorF vPos;
00202 StMuPrimaryVertex *muVert=mu->primaryVertex();
00203 if(muVert){
00204 vPos=muVert->position();
00205 }
00206 else if(mAUAU200){
00207 vPos.setX(muEvent->primaryVertexPosition().x());
00208 vPos.setY(muEvent->primaryVertexPosition().y());
00209 vPos.setZ(muEvent->primaryVertexPosition().z());
00210 }
00211 else{
00212 vPos.setX(0);
00213 vPos.setY(0);
00214 vPos.setZ(0);
00215 }
00216
00217
00218 if(mAUAU200){
00219 if(fabs(vPos.z())>30.){
00220 if(mDebug) cout<<"vertex out of range"<<endl;
00221 return kStOK;
00222 }
00223 }
00224
00225
00226 mTrigger=0;
00227 if(muEvent->triggerIdCollection().nominal().isTrigger(mTrig[0])||
00228 muEvent->triggerIdCollection().nominal().isTrigger(mTrig[1])) mTrigger+=1;
00229 if(muEvent->triggerIdCollection().nominal().isTrigger(mTrig[2])) mTrigger+=2;
00230 if(muEvent->triggerIdCollection().nominal().isTrigger(mTrig[3])) mTrigger+=4;
00231 if(muEvent->triggerIdCollection().nominal().isTrigger(2002)) mTrigger+=8;
00232
00233 if(mTrigger==0){
00234 if(mDebug) cout<<"nothing more useless than an event with no trigger!"<<endl;
00235 return kStOK;
00236 }
00237
00238 StEmcCollection *emcCol=(StEmcCollection*)mu->emcCollection();
00239 if(!emcCol){
00240 if(mDebug) cout<<"no StEmcCollection"<<endl;
00241 h_EvSum->Fill(8);
00242 return kStOK;
00243 }
00244 StEmcDetector *emcDet=(StEmcDetector*)emcCol->detector(kBarrelEmcTowerId);
00245 if(!emcDet)
00246 {
00247 if(mDebug) cout<<"no emcDet!!"<<endl;
00248 h_EvSum->Fill(9);
00249 return kStOK;
00250 }
00251
00252 Int_t hiTowId=-1;
00253 Int_t hiTowAdc=-1;
00254 Int_t hiTowStatus=-1;
00255 Float_t hiTowEnergy=0.;
00256 if(mTrigger>1){
00257 for(UInt_t m=1;m<=emcDet->numberOfModules();m++){
00258 StEmcModule *emcMod=emcDet->module(m);
00259 if(!emcMod){
00260 if(mDebug) cout<<"no emcMod for module "<<m<<endl;
00261 continue;
00262 }
00263 StSPtrVecEmcRawHit& modHits=emcMod->hits();
00264 for(UInt_t h=0;h<modHits.size();h++){
00265 Int_t adc=modHits[h]->adc();
00266 Float_t nrg=modHits[h]->energy();
00267 if(adc>hiTowAdc && nrg>0.){
00268 hiTowAdc=adc;
00269 hiTowEnergy=nrg;
00270 bemcGeom->getId(m,modHits[h]->eta(),modHits[h]->sub(),hiTowId);
00271 }
00272 }
00273 }
00274
00275 if(hiTowId>0) mBemcTables->getStatus(BTOW,hiTowId,hiTowStatus);
00276 }
00277
00278 Int_t nPrimTracks=0;
00279 Int_t nGlobalTracks=0;
00280 Int_t nGoodPrimaries=0;
00281 Int_t nGoodPrimBarrel=0;
00282 Int_t nGoodGlobals=0;
00283 Float_t TpcPt=0.;
00284 Float_t TpcPtBarrelWest=0.;
00285
00286 StMuTrack *ptrack;
00287 if(!mAUAU200){
00288 nPrimTracks=mu->primaryTracks()->GetEntries();
00289 for(Int_t i=0;i<nPrimTracks;i++){
00290 ptrack=(StMuTrack*)mu->primaryTracks()->UncheckedAt(i);
00291 if(ptrack->bad()) continue;
00292 if(ptrack->dca().mag()<3. && ptrack->nHitsFit()>15){
00293 nGoodPrimaries++;
00294 TpcPt+=ptrack->momentum().mag();
00295 if(fabs(ptrack->eta() - 0.5)<0.5){
00296 nGoodPrimBarrel++;
00297 TpcPtBarrelWest+=ptrack->momentum().mag();
00298 }
00299 }
00300 }
00301 StMuTrack *gtrack;
00302 nGlobalTracks=mu->globalTracks()->GetEntries();
00303 for(Int_t i=0;i<nGlobalTracks;i++){
00304 gtrack=(StMuTrack*)mu->globalTracks()->UncheckedAt(i);
00305 if(gtrack->bad()) continue;
00306 if(gtrack->dca().mag()<3. && gtrack->nHitsFit()>15)
00307 nGoodGlobals++;
00308 }
00309 }
00310
00311 Int_t ftpcRefMultTracks=muEvent->refMultFtpcEast();
00312 if(mAUAU200) ftpcRefMultTracks=muEvent->refMult();
00313 if(mDebug) cout<<"event has refMult: "<<ftpcRefMultTracks<<endl;
00314
00315
00316
00317 Float_t EnergyNeutral=0.;
00318
00319 StEmcPoint *pt=0;
00320 StSPtrVecEmcPoint& barrelPoints=emcCol->barrelPoints();
00321 for(UInt_t i=0;i<barrelPoints.size();i++){
00322
00323 pt=(StEmcPoint*)barrelPoints[i];
00324 Float_t dist=9999.;
00325 if(!calcDistanceTrackToPoint(pt,mu,dist)){
00326 cout<<"wrong bfield or no event pointer (shouldn't see this)"<<endl;
00327 }
00328 MyPoint *mypt=new MyPoint();
00329
00330 if(pt->position().pseudoRapidity()>0.0)
00331 EnergyNeutral+=pt->energy();
00332
00333 StPtrVecEmcCluster& etaClus=pt->cluster(kBarrelSmdEtaStripId);
00334 StPtrVecEmcCluster& phiClus=pt->cluster(kBarrelSmdPhiStripId);
00335 Int_t nEtaHits=0;
00336 Int_t nPhiHits=0;
00337 Float_t etaEnergy=0.;
00338 Float_t phiEnergy=0.;
00339 Float_t etaWidth=0.;
00340 Float_t phiWidth=0.;
00341 if(etaClus.size()>0){
00342 StEmcCluster *clE=(StEmcCluster*)etaClus[0];
00343 StPtrVecEmcRawHit& hE=clE->hit();
00344 nEtaHits=hE.size();
00345 etaEnergy=clE->energy();
00346 etaWidth=clE->sigmaEta();
00347 }
00348 if(phiClus.size()>0){
00349 StEmcCluster *clP=(StEmcCluster*)phiClus[0];
00350 StPtrVecEmcRawHit& hP=clP->hit();
00351 nPhiHits=hP.size();
00352 phiEnergy=clP->energy();
00353 phiWidth=clP->sigmaPhi();
00354 }
00355
00356 if( mAUAU200 && (nEtaHits<1||nPhiHits<1) ){
00357 delete mypt;
00358 continue;
00359 }
00360
00361 mypt->setEnergy(pt->energy());
00362 mypt->setQuality((Int_t)pt->quality());
00363 mypt->setPosition(pt->position().x(),pt->position().y(),pt->position().z());
00364 mypt->setDistanceToTrack(dist);
00365 mypt->setHitsEta(nEtaHits);
00366 mypt->setWidthEta(etaWidth);
00367 mypt->setEnergyEta(etaEnergy);
00368 mypt->setHitsPhi(nPhiHits);
00369 mypt->setWidthPhi(phiWidth);
00370 mypt->setEnergyPhi(phiEnergy);
00371
00372 mEvent->addPoint(mypt);
00373
00374 delete mypt;
00375 }
00376
00377
00378 mEvent->setPrescale(0,mPs_mb);
00379 mEvent->setPrescale(1,mPs_mb2);
00380 mEvent->setPrescale(2,mPs_ht1);
00381 mEvent->setPrescale(3,mPs_ht2);
00382 mEvent->setHighTowerAdc(hiTowAdc);
00383 mEvent->setHighTowerId(hiTowId);
00384 mEvent->setHighTowerStatus(hiTowStatus);
00385 mEvent->setHighTowerEnergy(hiTowEnergy);
00386 mEvent->setVertex(vPos.x(),vPos.y(),vPos.z());
00387 mEvent->setTrigger(mTrigger);
00388 mEvent->setGoodPrimaries(nGoodPrimaries);
00389 mEvent->setGoodPrimBarrel(nGoodPrimBarrel);
00390 mEvent->setGoodGlobals(nGoodGlobals);
00391 mEvent->setRefMult(ftpcRefMultTracks);
00392 mEvent->setEnergyInBarrel(EnergyNeutral);
00393 mEvent->setMomentumInTpc(TpcPt);
00394 mEvent->setMomentumInTpcWest(TpcPtBarrelWest);
00395 mEvent->setCtbSum(mips);
00396 mEvent->setBbcSumEast(bbcEast);
00397 mEvent->setBbcSumWest(bbcWest);
00398 mEvent->setZdcSumEast(zdcEast);
00399 mEvent->setZdcSumWest(zdcWest);
00400 mEvent->setZdcVertexZ(zdcVertex);
00401 mEvent->setBbcVertexZ(BBCVertexZ);
00402
00403 mEventTree->Fill();
00404
00405 return kStOK;
00406 }
00407
00408 Int_t StMyEventMaker::Finish()
00409 {
00410 saveHistograms();
00411 return kStOK;
00412 }
00413
00414 void StMyEventMaker::saveHistograms()
00415 {
00416 if(mDebug) cout<<"************ saving histograms ****************"<<endl;
00417 TFile *hfile=(TFile*)gROOT->FindObject(mFileName);
00418 if(hfile) hfile->Close();
00419 hfile=new TFile(mFileName,"RECREATE");
00420
00421 h_EvSum->Write();
00422 mEventTree->Write();
00423
00424 hfile->Close();
00425 }
00426
00427 Bool_t StMyEventMaker::calcDistanceTrackToPoint(StEmcPoint *point,StMuDst *currMuDst,Float_t &distanceToTrack)
00428 {
00429 StMuEvent *muEv=(StMuEvent*)currMuDst->event();
00430 if(!muEv) return kFALSE;
00431 Double_t bFld=muEv->magneticField()/10.;
00432 if(mDebug) cout<<"summary()->magneticField()="<<bFld<<"(Tesla)"<<endl;
00433
00434 if(fabs(bFld)<0.01){
00435 if(mDebug) cout<<"wrong/no BField!"<<endl;
00436 return kFALSE;
00437 }
00438
00439 StEmcPosition* emcPosition = new StEmcPosition();
00440 StThreeVectorD pos,mom;
00441 distanceToTrack=999.;
00442 StMuTrack *track;
00443
00444 Int_t ntracks=currMuDst->primaryTracks()->GetEntries();
00445 for(Int_t i=0;i<ntracks;i++){
00446 track=(StMuTrack*)currMuDst->primaryTracks()->UncheckedAt(i);
00447 if(track->bad()) continue;
00448
00449 StPhysicalHelixD helix = track->outerHelix();
00450 Bool_t projOK=emcPosition->projTrack(&pos,&mom,&helix,bFld,230.705,1);
00451 if(!projOK) continue;
00452 Float_t d=(pos - point->position()).mag();
00453
00454 if(d<distanceToTrack&&d>0.) distanceToTrack=d;
00455 }
00456
00457
00458 delete emcPosition;
00459 return kTRUE;
00460 }
00461
00462
00463
00464
00465
00466
