00001 #include <TROOT.h>
00002 #include <TTree.h>
00003 #include <TH2F.h>
00004 #include <TVector3.h>
00005 #include <TFile.h>
00006 #include <TMath.h>
00007 #include <assert.h>
00008
00009 #include <StEvent.h>
00010 #include <StEventTypes.h>
00011 #include <StEmcUtil/geometry/StEmcGeom.h>
00012 #include <StEmcUtil/projection/StEmcPosition.h>
00013 #include <StEmcUtil/database/StBemcTables.h>
00014 #include <StEvent/StEnumerations.h>
00015 #include <StEmcRawMaker/defines.h>
00016 #include <StEmcADCtoEMaker/StEmcADCtoEMaker.h>
00017
00018 #include <StMcEvent/StMcEvent.hh>
00019 #include <StMcEvent/StMcVertex.hh>
00020 #include <StMcEvent/StMcTrack.hh>
00021
00022
00023 #include <St_DataSet.h>
00024 #include <St_DataSetIter.h>
00025 #include <tables/St_g2t_event_Table.h>
00026 #include <tables/St_particle_Table.h>
00027 #include <tables/St_g2t_pythia_Table.h>
00028
00029 #include <StDetectorDbMaker/StDetectorDbTriggerID.h>
00030 #include <St_db_Maker/St_db_Maker.h>
00031
00032 #include <StEmcTriggerMaker/StEmcTriggerMaker.h>
00033
00034
00035 #include <StEmcPool/StPhotonCommon/MyEvent.h>
00036 #include <StEmcPool/StPhotonCommon/MyPoint.h>
00037 #include <StEmcPool/StPhotonCommon/MyMcTrack.h>
00038
00039 #include "StPhotonMaker.h"
00040
00041 ClassImp(StPhotonMaker)
00042
00043 StPhotonMaker::StPhotonMaker(const char *name,const char *outputFile,
00044 const char *type,const char *coll,Bool_t debug):StMaker(name)
00045 {
00046 mFileName=outputFile;
00047
00048 mMc=kFALSE;
00049 mEmbed=kFALSE;
00050 mReal=kFALSE;
00051 mPythia=kFALSE;
00052 mHijing=kFALSE;
00053 if(strcmp(type,"mc")==0) mMc=kTRUE;
00054 else if(strcmp(type,"embed")==0) mEmbed=kTRUE;
00055 else if(strcmp(type,"pythia")==0) mPythia=kTRUE;
00056 else if(strcmp(type,"hijing")==0) mHijing=kTRUE;
00057 else mReal=kTRUE;
00058
00059 mDAU=kFALSE;
00060 mPP04=kFALSE;
00061 mPP05=kFALSE;
00062 if(strcmp(coll,"dAu")==0){
00063 mDAU=kTRUE;
00064 mTrig[0]=2001;
00065 mTrig[1]=2003;
00066 mTrig[2]=2201;
00067 mTrig[3]=2202;
00068 }
00069 if(strcmp(coll,"pp04")==0){
00070 mPP04=kTRUE;
00071 mTrig[0]=45010;
00072 mTrig[1]=10;
00073 mTrig[2]=45201;
00074 mTrig[3]=45202;
00075 }
00076 if(strcmp(coll,"pp05")==0){
00077 mPP05=kTRUE;
00078 mTrig[0]=96011;
00079 mTrig[1]=-10000;
00080 mTrig[2]=96201;
00081 mTrig[3]=96211;
00082 }
00083 mDebug=debug;
00084
00085 mBemcTables=new StBemcTables();
00086 mEvent=new MyEvent();
00087 }
00088
00089 StPhotonMaker::~StPhotonMaker()
00090 {
00091
00092 }
00093
00094 Int_t StPhotonMaker::Init()
00095 {
00096
00097 h_EvSum=new TH1F("h_EvSum","see code for details",22,-1.5,20.5);
00098 h_bsmdeAdc=new TH1F("h_bsmdeAdc","adc",1200,-0.5,1199.5);
00099 h_bsmdpAdc=new TH1F("h_bsmdpAdc","adc",1200,-0.5,1199.5);
00100 h_btowAdc=new TH1F("h_btowAdc","adc",1200,-0.5,1199.5);
00101 h_bsmdeEn=new TH1F("h_bsmdeEn","energy",55,-1.,10.);
00102 h_btowEn=new TH1F("h_btowEn","energy",55,-1.,10.);
00103 h_btowEnVsAdc=new TH2F("h_btowEnVsAdc","h_btowEnVsAdc",1200,-0.5,1199.5,100,0.,30.);
00104
00105 mEventTree=new TTree("mEventTree","tree with my event structure");
00106 mEventTree->Branch("branch","MyEvent",&mEvent);
00107
00108 mN=0;
00109 mRunId=0;
00110 mRunPrev=0;
00111 mEventId=0;
00112 mTrigger=0;
00113 mDate=0;
00114 mTime=0;
00115 mPs_mb=0;
00116 mPs_mb2=0;
00117 mPs_ht1=0;
00118 mPs_ht2=0;
00119
00120 return StMaker::Init();
00121 }
00122
00123 Int_t StPhotonMaker::Make()
00124 {
00125 mN++;
00126 if(mDebug) cout<<endl<<"#### processing event "<<mN<<endl<<endl;
00127 h_EvSum->Fill(1);
00128
00129 StEmcGeom *bemcGeom=StEmcGeom::getEmcGeom("bemc");
00130
00131 StEvent *event=0;
00132 event = (StEvent*)GetInputDS("StEvent");
00133 if(!event){
00134 h_EvSum->Fill(2);
00135 if(mDebug) cout << "************ StMyPhotonMaker::Make() no event pointer!! *********" << endl;
00136 return kStOK;
00137 }
00138
00139
00140
00141
00142
00143
00144
00145 if(!mHijing && !mMc && !mPythia && mAdcMaker->isCorrupted()){
00146 if(mDebug) cout<<"event corrupted!"<<endl;
00147 h_EvSum->Fill(-1);
00148 return kStOK;
00149 }
00150
00151 mRunPrev=mRunId;
00152 if(!mMc){
00153 mRunId=event->runId();
00154 mEventId=event->id();
00155 }
00156 else{
00157 mEventId=mN;
00158 mRunId=0;
00159 }
00160 if(mDebug) cout<<"starting run: "<<mRunId<<" and event: "<<mEventId<<endl;
00161
00162
00163 if(mRunId!=mRunPrev&&mReal){
00164 StDetectorDbTriggerID& DetDbTrigId = *(StDetectorDbTriggerID::instance());
00165 for(UInt_t i=0;i<DetDbTrigId.getL0NumRows();i++){
00166 if(mDebug) cout<<"prescales: "<<DetDbTrigId.getPsL0(i)<<endl<<endl;
00167 if(DetDbTrigId.getL0OfflineTrgId(i)==mTrig[0]) mPs_mb=DetDbTrigId.getPsL0(i);
00168 if(DetDbTrigId.getL0OfflineTrgId(i)==mTrig[1]) mPs_mb2=DetDbTrigId.getPsL0(i);
00169 if(DetDbTrigId.getL0OfflineTrgId(i)==mTrig[2]) mPs_ht1=DetDbTrigId.getPsL0(i);
00170 if(DetDbTrigId.getL0OfflineTrgId(i)==mTrig[3]) mPs_ht2=DetDbTrigId.getPsL0(i);
00171 }
00172 mBemcTables->loadTables((StMaker*)this);
00173 }
00174
00175 mDate=(Int_t)mDbMaker->GetDateTime().GetDate();
00176 mTime=(Int_t)mDbMaker->GetDateTime().GetTime();
00177
00178
00179 Float_t mips=-1.;
00180 Float_t zdcEast=0.;
00181 Float_t zdcWest=0.;
00182 Float_t zdcVertex=0.;
00183 Float_t bbcEast=0.;
00184 Float_t bbcWest=0.;
00185 StTriggerDetectorCollection *trigDetColl=event->triggerDetectorCollection();
00186 if(trigDetColl){
00187 StCtbTriggerDetector &CTB=trigDetColl->ctb();
00188 mips=CTB.mips(0);
00189 StZdcTriggerDetector &ZDC=trigDetColl->zdc();
00190 zdcEast=ZDC.adcSum(east);
00191 zdcWest=ZDC.adcSum(west);
00192 zdcVertex=ZDC.vertexZ();
00193 StBbcTriggerDetector &BBC=trigDetColl->bbc();
00194 bbcEast=BBC.adcSumEast();
00195 bbcWest=BBC.adcSumWest();
00196 }
00197
00198
00199 mEvent=new MyEvent(mRunId,mEventId,mDate,mTime);
00200 if(!mEvent){
00201 if(mDebug) cout<<"couldn't create MyEvent* !!"<<endl;
00202 return kStErr;
00203 }
00204
00205
00206 StPrimaryVertex *pVert=event->primaryVertex();
00207 StMcVertex *mc_pVert=0;
00208 StMcEvent* mc_event=0;
00209 if(mPythia||mHijing){
00210 mc_event=(StMcEvent*)GetInputDS("StMcEvent");
00211 if(!mc_event){
00212 h_EvSum->Fill(4);
00213 if(mDebug) cout<<"StPhotonMaker::Make() no McEvent"<<endl;
00214 return kStOK;
00215 }
00216 mc_pVert=mc_event->primaryVertex();
00217 if(!mc_pVert){
00218 h_EvSum->Fill(5);
00219 if(mDebug) cout<<"StPhotonMaker::Make() no mc vertex!"<<endl;
00220 return kStOK;
00221 }
00222 }
00223 if(mMc||mEmbed){
00224 mc_event=(StMcEvent*)GetInputDS("StMcEvent");
00225 if(!mc_event){
00226 h_EvSum->Fill(4);
00227 if(mDebug) cout<<"StPhotonMaker::Make() no McEvent"<<endl;
00228 return kStOK;
00229 }
00230 mc_pVert=mc_event->primaryVertex();
00231 if(!mc_pVert){
00232 h_EvSum->Fill(5);
00233 if(mDebug) cout<<"StPhotonMaker::Make() no mc vertex!"<<endl;
00234 return kStOK;
00235 }
00236
00237 StPtrVecMcTrack& trcks=mc_pVert->daughters();
00238 Int_t nDaughters=mc_pVert->numberOfDaughters();
00239 if(mDebug) cout<<"number of MC part. from vertex: "<<nDaughters<<endl;
00240 if(nDaughters==1){
00241 StMcTrackIterator it=trcks.begin();
00242 Int_t pid=(*it)->particleDefinition()->pdgEncoding();
00243 if(mDebug) cout<<"daughter pid: "<<pid<<endl;
00244
00245
00246 if(pid==22){
00247 MyMcTrack *mytr=new MyMcTrack();
00248 mytr->setId(pid);
00249 mytr->setEnergy((*it)->energy());
00250 if((*it)->stopVertex()){
00251 if((*it)->stopVertex()->daughters().size()){
00252 mytr->setStopRadius((*it)->stopVertex()->position().perp());
00253 if(mDebug) cout<<"radius: "<<mytr->stopRadius()<<endl;
00254 }
00255 else if(mDebug) cout<<"no size"<<endl;
00256 }
00257 else if(mDebug) cout<<"no stop"<<endl;
00258 mytr->setMomentum((*it)->momentum().x(),(*it)->momentum().y(),(*it)->momentum().z());
00259
00260 StThreeVectorF v;
00261 v.setX(bemcGeom->Radius()*cos((*it)->momentum().phi()));
00262 v.setY(bemcGeom->Radius()*sin((*it)->momentum().phi()));
00263 v.setZ(bemcGeom->Radius()/tan((*it)->momentum().theta()));
00264 v=v + mc_pVert->position();
00265
00266 mytr->setPosition(v.x(),v.y(),v.z());
00267
00268 mEvent->setMcTrack(mytr);
00269 delete mytr;
00270 }
00271 else if(pid==111||pid==221||pid==2112||pid==-2112||pid==211||pid==130){
00272
00273 MyMcTrack *mytr=new MyMcTrack();
00274 mytr->setId(pid);
00275 mytr->setEnergy((*it)->energy());
00276
00277 if((*it)->stopVertex()){
00278 if((*it)->stopVertex()->daughters().size()){
00279 mytr->setStopRadius((*it)->stopVertex()->position().perp());
00280 }
00281 }
00282
00283 StThreeVectorF v;
00284 v.setX(bemcGeom->Radius()*cos((*it)->momentum().phi()));
00285 v.setY(bemcGeom->Radius()*sin((*it)->momentum().phi()));
00286 v.setZ(bemcGeom->Radius()/tan((*it)->momentum().theta()));
00287 v=v + mc_pVert->position();
00288
00289 mytr->setMomentum((*it)->momentum().x(),(*it)->momentum().y(),(*it)->momentum().z());
00290 mytr->setPosition(v.x(),v.y(),v.z());
00291
00292 mEvent->setMcTrack(mytr);
00293 delete mytr;
00294
00295
00296 if((*it)->stopVertex()){
00297 StMcVertex *sVert=(*it)->stopVertex();
00298 StPtrVecMcTrack& gammaTracks=sVert->daughters();
00299 for(StMcTrackIterator itd=gammaTracks.begin();itd!=gammaTracks.end();itd++){
00300 Int_t pid2=(*itd)->particleDefinition()->pdgEncoding();
00301 MyMcTrack *mytr2=new MyMcTrack();
00302 mytr2->setId(pid2);
00303 mytr2->setEnergy((*itd)->energy());
00304
00305 if((*itd)->stopVertex()){
00306 if((*itd)->stopVertex()->daughters().size()){
00307 mytr2->setStopRadius((*it)->stopVertex()->position().perp());
00308 }
00309 }
00310
00311 StThreeVectorF vv;
00312 vv.setX(bemcGeom->Radius()*cos((*itd)->momentum().phi()));
00313 vv.setY(bemcGeom->Radius()*sin((*itd)->momentum().phi()));
00314 vv.setZ(bemcGeom->Radius()/tan((*itd)->momentum().theta()));
00315 vv=vv + mc_pVert->position();
00316
00317 mytr2->setMomentum((*itd)->momentum().x(),(*itd)->momentum().y(),(*itd)->momentum().z());
00318 mytr2->setPosition(vv.x(),vv.y(),vv.z());
00319
00320 mEvent->addMcPhoton(mytr2);
00321 delete mytr2;
00322 }
00323 }
00324 }
00325 else
00326 cout<<"what kind of MC particle is this??? (shouldn't see this)"<<endl;
00327 }
00328 else if(mMc && nDaughters>1){
00329 if(mDebug) cout<<"more than one particle per event"<<endl;
00330
00331 StMcTrackIterator it=trcks.begin();
00332 Int_t pid=(*it)->particleDefinition()->pdgEncoding();
00333 if(mDebug) cout<<"daughter pid: "<<pid<<endl;
00334 if(pid==-2112 || pid==2112 || pid==221){
00335 for(UInt_t i=0;i<trcks.size();i++){
00336 Float_t neutronenergy=trcks[i]->energy();
00337 MyMcTrack *myTR=new MyMcTrack();
00338 myTR->setId(pid);
00339 myTR->setEnergy(neutronenergy);
00340 myTR->setMomentum(trcks[i]->momentum().x(),trcks[i]->momentum().y(),trcks[i]->momentum().z());
00341 StThreeVectorF v;
00342 v.setX(bemcGeom->Radius()*cos(trcks[i]->momentum().phi()));
00343 v.setY(bemcGeom->Radius()*sin(trcks[i]->momentum().phi()));
00344 v.setZ(bemcGeom->Radius()/tan(trcks[i]->momentum().theta()));
00345 v=v + mc_pVert->position();
00346 myTR->setPosition(v.x(),v.y(),v.z());
00347 mEvent->addMcPion(myTR);
00348 delete myTR;
00349 }
00350 }
00351 }
00352 }
00353
00354 if(!pVert&&!mc_pVert){
00355 h_EvSum->Fill(6);
00356 if(mDebug) cout<<"StPhotonMaker::Make() no primary vertex"<<endl;
00357 return kStOK;
00358 }
00359
00360 if(mPythia){
00361
00362 TDataSet *gEvent=GetDataSet("geant");
00363 TDataSetIter geantDstI(gEvent);
00364
00365
00366
00367
00368
00369
00370 St_g2t_pythia *Pg2t_pythia=(St_g2t_pythia*)geantDstI("g2t_pythia");
00371 if(Pg2t_pythia){
00372 g2t_pythia_st *g2t_pythia1=Pg2t_pythia->GetTable();
00373 if(g2t_pythia1){
00374 float hard_p=g2t_pythia1->hard_p;
00375 mEvent->setPartonPt(hard_p);
00376 }
00377 }
00378 }
00379 if(mPythia||mHijing){
00380 StPtrVecMcTrack& tracs=mc_event->tracks();
00381 for(UInt_t i=0;i<tracs.size();i++){
00382
00383 if(tracs[i]->geantId()==7){
00384 Float_t pionenergy=tracs[i]->energy();
00385 MyMcTrack *myTR=new MyMcTrack();
00386 myTR->setId(111);
00387 myTR->setEnergy(pionenergy);
00388 myTR->setMomentum(tracs[i]->momentum().x(),tracs[i]->momentum().y(),tracs[i]->momentum().z());
00389 StThreeVectorF v;
00390 v.setX(bemcGeom->Radius()*cos(tracs[i]->momentum().phi()));
00391 v.setY(bemcGeom->Radius()*sin(tracs[i]->momentum().phi()));
00392 v.setZ(bemcGeom->Radius()/tan(tracs[i]->momentum().theta()));
00393 v=v + mc_pVert->position();
00394 myTR->setPosition(v.x(),v.y(),v.z());
00395 mEvent->addMcPion(myTR);
00396 delete myTR;
00397 }
00398 if(tracs[i]->geantId()==1){
00399 Float_t photonenergy=tracs[i]->energy();
00400 MyMcTrack *myTR=new MyMcTrack();
00401 myTR->setId(22);
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412 myTR->setEnergy(photonenergy);
00413 myTR->setMomentum(tracs[i]->momentum().x(),tracs[i]->momentum().y(),tracs[i]->momentum().z());
00414 StThreeVectorF v;
00415 v.setX(bemcGeom->Radius()*cos(tracs[i]->momentum().phi()));
00416 v.setY(bemcGeom->Radius()*sin(tracs[i]->momentum().phi()));
00417 v.setZ(bemcGeom->Radius()/tan(tracs[i]->momentum().theta()));
00418 v=v + mc_pVert->position();
00419 myTR->setPosition(v.x(),v.y(),v.z());
00420 mEvent->addMcPhoton(myTR);
00421 delete myTR;
00422 }
00423 }
00424 }
00425
00426 if(!pVert&&mEmbed){
00427 h_EvSum->Fill(7);
00428 if(mDebug) cout<<"StPhotonMaker::Make() no primary vertex ****"<<endl;
00429 return kStOK;
00430 }
00431 StThreeVectorF vPos,mc_vPos;
00432 if(mMc||mEmbed||mPythia||mHijing){
00433 mc_vPos=mc_pVert->position();
00434 if(mMc||mPythia) vPos=mc_vPos;
00435 }
00436 if(mEmbed||mReal||mHijing){
00437 if(pVert) vPos=pVert->position();
00438 }
00439
00440
00441 mTrigger=0;
00442 if(mReal&&event->triggerIdCollection()&&event->triggerIdCollection()->nominal()){
00443 if(event->triggerIdCollection()->nominal()->isTrigger(mTrig[0])||
00444 event->triggerIdCollection()->nominal()->isTrigger(mTrig[1])) mTrigger+=1;
00445 if(event->triggerIdCollection()->nominal()->isTrigger(mTrig[2])) mTrigger+=2;
00446 if(event->triggerIdCollection()->nominal()->isTrigger(mTrig[3])) mTrigger+=4;
00447 if(event->triggerIdCollection()->nominal()->isTrigger(2002)) mTrigger+=8;
00448 }
00449
00450 StEmcTriggerMaker *triggermaker=dynamic_cast<StEmcTriggerMaker*>(GetMaker("bemctrigger"));
00451 assert(triggermaker);
00452
00453 if(mMc||mEmbed){
00454 mTrigger=1;
00455 if(mDate<20060001){
00456 if(triggermaker->is2005HT1()) mTrigger+=2;
00457 if(triggermaker->is2005HT2()) mTrigger+=4;
00458 }
00459 else{
00460 cout<<"check timestamp -> triggermaker"<<endl;
00461 assert(0);
00462 }
00463 }
00464 if(mPythia){
00465 mTrigger=1;
00466
00467 if(mDate<20050101) assert(0);
00468 if(triggermaker->is2005HT1()) mTrigger+=2;
00469 if(triggermaker->is2005HT2()) mTrigger+=4;
00470 }
00471 if(mHijing){
00472 mTrigger=1;
00473 }
00474
00475 if(mTrigger==0){
00476 if(mDebug) cout<<"nothing more useless than an event with no trigger!"<<endl;
00477 return kStOK;
00478 }
00479
00480 StEmcCollection *emcCol=(StEmcCollection*)event->emcCollection();
00481 if(!emcCol){
00482 if(mDebug) cout<<"no emccollection!!"<<endl;
00483 h_EvSum->Fill(8);
00484 return kStOK;
00485 }
00486 StEmcDetector *emcDet=(StEmcDetector*)emcCol->detector(kBarrelEmcTowerId);
00487 if(!emcDet){
00488 if(mDebug) cout<<"no emcDet!!"<<endl;
00489 h_EvSum->Fill(9);
00490 return kStOK;
00491 }
00492
00493 Int_t hiTowId=-1;
00494 Int_t hiTowAdc=-1;
00495 Int_t hiTowStatus=-1;
00496 Float_t hiTowEnergy=0.;
00497 if(mTrigger>1 || mEmbed || mMc){
00498 for(UInt_t m=1;m<=emcDet->numberOfModules();m++){
00499 StEmcModule *emcMod=emcDet->module(m);
00500 if(!emcMod){
00501 if(mDebug) cout<<"no emcMod for module "<<m<<endl;
00502 continue;
00503 }
00504 StSPtrVecEmcRawHit& modHits=emcMod->hits();
00505 for(UInt_t h=0;h<modHits.size();h++){
00506 Int_t adc=modHits[h]->adc();
00507 Float_t nrg=modHits[h]->energy();
00508
00509
00510
00511 if(adc>hiTowAdc && nrg>0.){
00512 hiTowAdc=adc;
00513 hiTowEnergy=nrg;
00514 bemcGeom->getId(m,modHits[h]->eta(),modHits[h]->sub(),hiTowId);
00515 }
00516 }
00517 }
00518
00519 if(hiTowId>0) mBemcTables->getStatus(BTOW,hiTowId,hiTowStatus);
00520 }
00521
00522 StEmcDetector *bsmdeDet=(StEmcDetector*)emcCol->detector(kBarrelSmdEtaStripId);
00523 if(!bsmdeDet){
00524 return kStOK;
00525 }
00526 StEmcDetector *bsmdpDet=(StEmcDetector*)emcCol->detector(kBarrelSmdPhiStripId);
00527 if(!bsmdpDet){
00528 return kStOK;
00529 }
00530 for(UInt_t m2=1;m2<=bsmdeDet->numberOfModules();m2++){
00531 StEmcModule *bsmdeMod=bsmdeDet->module(m2);
00532 StSPtrVecEmcRawHit& bsmdeHits=bsmdeMod->hits();
00533 for(UInt_t h2=0;h2<bsmdeHits.size();h2++){
00534 h_bsmdeAdc->Fill(bsmdeHits[h2]->adc());
00535 }
00536 }
00537 for(UInt_t m2=1;m2<=bsmdpDet->numberOfModules();m2++){
00538 StEmcModule *bsmdpMod=bsmdpDet->module(m2);
00539 StSPtrVecEmcRawHit& bsmdpHits=bsmdpMod->hits();
00540 for(UInt_t h2=0;h2<bsmdpHits.size();h2++){
00541 h_bsmdpAdc->Fill(bsmdpHits[h2]->adc());
00542 }
00543 }
00544
00545
00546 StSPtrVecEmcPoint& barrelPoints=emcCol->barrelPoints();
00547 Int_t nPrimTracks=0;
00548 Int_t nGoodPrimaries=0;
00549 Int_t nGoodPrimBarrel=0;
00550 Int_t nGoodGlobals=0;
00551 Float_t TpcPt=0.;
00552 Float_t TpcPtBarrelWest=0.;
00553 if(event->primaryVertex()){
00554 nPrimTracks=pVert->numberOfDaughters();
00555 StSPtrVecTrackNode& trackNode=event->trackNodes();
00556 if(mDebug) cout<<"number of tracknodes: "<<trackNode.size()<<endl;
00557 for(UInt_t j=0;j<trackNode.size();j++){
00558 StPrimaryTrack* track=static_cast<StPrimaryTrack*>(trackNode[j]->track(primary));
00559 if(track && track->flag()>0 && track->geometry()){
00560
00561 Double_t pathlength=track->geometry()->helix().pathLength(vPos,kFALSE);
00562 StThreeVectorD v=track->geometry()->helix().at(pathlength) - vPos;
00563 Float_t dca=v.mag();
00564
00565 const StTrackFitTraits& fitTraits=track->fitTraits();
00566 unsigned short numFitPoints=fitTraits.numberOfFitPoints();
00567 if(numFitPoints>15&&dca<3.){
00568 nGoodPrimaries++;
00569 Float_t tracEta=track->geometry()->momentum().pseudoRapidity();
00570 if(tracEta>0.&&tracEta<1.){
00571 nGoodPrimBarrel++;
00572 TpcPtBarrelWest+=track->geometry()->momentum().perp();
00573 }
00574 TpcPt+=track->geometry()->momentum().perp();
00575 }
00576 }
00577 StGlobalTrack* gtrack=static_cast<StGlobalTrack*>(trackNode[j]->track(global));
00578 if(gtrack && gtrack->flag()>0 && gtrack->geometry()){
00579 const StTrackFitTraits& fitTraitsG=gtrack->fitTraits();
00580 unsigned short numFitPointsG=fitTraitsG.numberOfFitPoints();
00581 if(numFitPointsG>15) nGoodGlobals++;
00582 }
00583 }
00584 }
00585
00586
00587 Int_t ftpcRefMultTracks=0;
00588 if(event->primaryVertex()){
00589 const StSPtrVecPrimaryTrack& tracks=pVert->daughters();
00590 for(StSPtrVecPrimaryTrackConstIterator iter = tracks.begin(); iter != tracks.end(); iter++){
00591 StTrack* track = (*iter);
00592 if(!track->geometry()){
00593 if(mDebug) cout<<"no track geometry! -> NEXT"<<endl;
00594 continue;
00595 }
00596
00597 if(track->fitTraits().numberOfFitPoints()<6||(track->fitTraits().numberOfFitPoints()>11))
00598 continue;
00599
00600 if(track->geometry()->momentum().pseudoRapidity()>-2.8||track->geometry()->momentum().pseudoRapidity()<=-3.8)
00601 continue;
00602
00603 if(track->geometry()->momentum().perp()>=3.)
00604 continue;
00605
00606 StTrack *glt=track->node()->track(global);
00607 if(!glt){
00608 if(mDebug) cout<<"no global track!"<<endl;
00609 continue;
00610 }
00611 if(!glt->geometry()){
00612 if(mDebug) cout<<"no geom.!"<<endl;
00613 continue;
00614 }
00615 if(glt->geometry()->helix().distance(vPos)<3.) ftpcRefMultTracks++;
00616 }
00617 }
00618
00619
00620 StEmcPoint *pt=0;
00621 Float_t EnergyNeutral=0.;
00622 for(UInt_t i=0;i<barrelPoints.size();i++)
00623 {
00624 MyPoint *mypt=new MyPoint();
00625 pt=(StEmcPoint*)barrelPoints[i];
00626
00627 Float_t dist=9999.;
00628 if(!calcDistanceTrackToPoint(pt,dist)){
00629 cout<<"wrong bfield or no StEvent pointer (shouldn't see this)"<<endl;
00630 }
00631 if(pt->position().pseudoRapidity()>0.0)
00632 EnergyNeutral+=pt->energy();
00633
00634 StPtrVecEmcCluster& etaClus=pt->cluster(kBarrelSmdEtaStripId);
00635 StPtrVecEmcCluster& phiClus=pt->cluster(kBarrelSmdPhiStripId);
00636 Int_t nEtaHits=0;
00637 Int_t nPhiHits=0;
00638 Float_t etaEnergy=0.;
00639 Float_t phiEnergy=0.;
00640 Float_t etaWidth=0.;
00641 Float_t phiWidth=0.;
00642 if(etaClus.size()>0){
00643 StEmcCluster *clE=(StEmcCluster*)etaClus[0];
00644 StPtrVecEmcRawHit& hE=clE->hit();
00645 nEtaHits=hE.size();
00646 etaEnergy=clE->energy();
00647 etaWidth=clE->sigmaEta();
00648 }
00649 if(phiClus.size()>0){
00650 StEmcCluster *clP=(StEmcCluster*)phiClus[0];
00651 StPtrVecEmcRawHit& hP=clP->hit();
00652 nPhiHits=hP.size();
00653 phiEnergy=clP->energy();
00654 phiWidth=clP->sigmaPhi();
00655 }
00656 mypt->setEnergy(pt->energy());
00657 mypt->setQuality((Int_t)pt->quality());
00658 mypt->setPosition(pt->position().x(),pt->position().y(),pt->position().z());
00659 mypt->setDistanceToTrack(dist);
00660 mypt->setHitsEta(nEtaHits);
00661 mypt->setWidthEta(etaWidth);
00662 mypt->setEnergyEta(etaEnergy);
00663 mypt->setHitsPhi(nPhiHits);
00664 mypt->setWidthPhi(phiWidth);
00665 mypt->setEnergyPhi(phiEnergy);
00666 mEvent->addPoint(mypt);
00667
00668 delete mypt;
00669 }
00670
00671
00672
00673 mEvent->setPrescale(0,mPs_mb);
00674 mEvent->setPrescale(1,mPs_mb2);
00675 mEvent->setPrescale(2,mPs_ht1);
00676 mEvent->setPrescale(3,mPs_ht2);
00677 mEvent->setHighTowerAdc(hiTowAdc);
00678 mEvent->setHighTowerId(hiTowId);
00679 mEvent->setHighTowerStatus(hiTowStatus);
00680 mEvent->setHighTowerEnergy(hiTowEnergy);
00681 mEvent->setVertex(vPos.x(),vPos.y(),vPos.z());
00682 mEvent->setTrigger(mTrigger);
00683 mEvent->setGoodPrimaries(nGoodPrimaries);
00684 mEvent->setGoodPrimBarrel(nGoodPrimBarrel);
00685 mEvent->setGoodGlobals(nGoodGlobals);
00686 mEvent->setRefMult(ftpcRefMultTracks);
00687 mEvent->setEnergyInBarrel(EnergyNeutral);
00688 mEvent->setMomentumInTpc(TpcPt);
00689 mEvent->setMomentumInTpcWest(TpcPtBarrelWest);
00690 mEvent->setCtbSum(mips);
00691 mEvent->setBbcSumEast(bbcEast);
00692 mEvent->setBbcSumWest(bbcWest);
00693 mEvent->setZdcSumEast(zdcEast);
00694 mEvent->setZdcSumWest(zdcWest);
00695 mEvent->setZdcVertexZ(zdcVertex);
00696 if(mHijing){
00697 mEvent->setZdcVertexZ(mc_vPos.z());
00698 }
00699
00700 mEventTree->Fill();
00701
00702 return kStOK;
00703 }
00704
00705
00706
00707
00708
00709 Int_t StPhotonMaker::Finish()
00710 {
00711 saveHistograms();
00712 return kStOK;
00713 }
00714
00715
00716 void StPhotonMaker::saveHistograms()
00717 {
00718 if(mDebug) cout<<"************ saving histograms ****************"<<endl;
00719 TFile *hfile=(TFile*)gROOT->FindObject(mFileName);
00720 if(hfile) hfile->Close();
00721 hfile=new TFile(mFileName,"RECREATE");
00722
00723
00724 h_EvSum->Write();
00725 h_bsmdeAdc->Write();
00726 h_bsmdpAdc->Write();
00727
00728
00729
00730
00731
00732 mEventTree->Write();
00733
00734
00735 hfile->Close();
00736 }
00737
00738 void StPhotonMaker::setDbMaker(St_db_Maker *dbM)
00739 {
00740 mDbMaker=dbM;
00741 }
00742 void StPhotonMaker::setAdcMaker(StEmcADCtoEMaker *adc)
00743 {
00744 mAdcMaker=adc;
00745 }
00746 Bool_t StPhotonMaker::calcDistanceTrackToPoint(StEmcPoint *point,Float_t &distanceToTrack)
00747 {
00748 StEvent *event = (StEvent*)this->GetInputDS("StEvent");
00749 if(!event){
00750 if(mDebug) cout<<"can't get Event pointer"<<endl;
00751 return kFALSE;
00752 }
00753
00754 Double_t bFld=0.;
00755 StEventSummary* summary = event->summary();
00756 if(summary){
00757 bFld = summary->magneticField()/10.;
00758 if(mDebug) cout<<"summary()->magneticField()="<<bFld<<"(Tesla)"<<endl;
00759 }
00760 if(TMath::Abs(bFld)<0.01&&!mMc){
00761 if(mDebug) cout<<"wrong mBField!"<<endl;
00762 return kFALSE;
00763 }
00764 StEmcPosition* emcPosition = new StEmcPosition();
00765 StThreeVectorD pos, mom;
00766 StSPtrVecTrackNode& trackNodes = event->trackNodes();
00767 StTrack* track;
00768 distanceToTrack=999.;
00769 for (size_t nodeIndex=0; nodeIndex<trackNodes.size();nodeIndex++){
00770 size_t numberOfTracksInNode=trackNodes[nodeIndex]->entries(global);
00771 for(size_t trackIndex=0; trackIndex<numberOfTracksInNode;trackIndex++){
00772 track = trackNodes[nodeIndex]->track(primary,trackIndex);
00773 if (track && track->flag()>=0){
00774 if(track->geometry()->momentum().mag()<.5) continue;
00775 StPhysicalHelixD helix=track->outerGeometry()->helix();
00776 Bool_t projOK=emcPosition->projTrack(&pos,&mom,track,bFld,230.705,1);
00777 if(!projOK) continue;
00778 Float_t d=(pos - point->position()).mag();
00779 if(d<distanceToTrack&&d>0.) distanceToTrack=d;
00780 }
00781 }
00782 }
00783 delete emcPosition;
00784 return kTRUE;
00785 }
