StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StPhotonMaker.cxx
1 #include <TROOT.h>
2 #include <TTree.h>
3 #include <TH2F.h>
4 #include <TVector3.h>
5 #include <TFile.h>
6 #include <TMath.h>
7 #include <assert.h>
8 
9 #include <StEvent.h>
10 #include <StEventTypes.h>
11 #include <StEmcUtil/geometry/StEmcGeom.h>
12 #include <StEmcUtil/projection/StEmcPosition.h>
13 #include <StEmcUtil/database/StBemcTables.h>
14 #include <StEvent/StEnumerations.h>
15 #include <StEmcRawMaker/defines.h>
16 #include <StEmcADCtoEMaker/StEmcADCtoEMaker.h>
17 
18 #include <StMcEvent/StMcEvent.hh>
19 #include <StMcEvent/StMcVertex.hh>
20 #include <StMcEvent/StMcTrack.hh>
21 
22 //pythia record
23 #include <St_DataSet.h>
24 #include <St_DataSetIter.h>
25 #include <tables/St_g2t_event_Table.h>
26 #include <tables/St_particle_Table.h>
27 #include <tables/St_g2t_pythia_Table.h>
28 
29 #include <StDetectorDbMaker/StDetectorDbTriggerID.h>
30 #include <St_db_Maker/St_db_Maker.h>
31 
32 #include <StEmcTriggerMaker/StEmcTriggerMaker.h>
33 
34 //MyEvent
35 #include <StEmcPool/StPhotonCommon/MyEvent.h>
36 #include <StEmcPool/StPhotonCommon/MyPoint.h>
37 #include <StEmcPool/StPhotonCommon/MyMcTrack.h>
38 
39 #include "StPhotonMaker.h"
40 
41 ClassImp(StPhotonMaker)
42 
43  StPhotonMaker::StPhotonMaker(const char *name,const char *outputFile,
44  const char *type,const char *coll,Bool_t debug):StMaker(name)
45 {
46  mFileName=outputFile;
47 
48  mMc=kFALSE;
49  mEmbed=kFALSE;
50  mReal=kFALSE;
51  mPythia=kFALSE;
52  mHijing=kFALSE;
53  if(strcmp(type,"mc")==0) mMc=kTRUE;
54  else if(strcmp(type,"embed")==0) mEmbed=kTRUE;
55  else if(strcmp(type,"pythia")==0) mPythia=kTRUE;
56  else if(strcmp(type,"hijing")==0) mHijing=kTRUE;
57  else mReal=kTRUE;
58 
59  mDAU=kFALSE;
60  mPP04=kFALSE;
61  mPP05=kFALSE;
62  if(strcmp(coll,"dAu")==0){
63  mDAU=kTRUE;
64  mTrig[0]=2001;
65  mTrig[1]=2003;
66  mTrig[2]=2201;
67  mTrig[3]=2202;
68  }
69  if(strcmp(coll,"pp04")==0){
70  mPP04=kTRUE;
71  mTrig[0]=45010;
72  mTrig[1]=10;
73  mTrig[2]=45201;
74  mTrig[3]=45202;
75  }
76  if(strcmp(coll,"pp05")==0){
77  mPP05=kTRUE;
78  mTrig[0]=96011;
79  mTrig[1]=-10000;
80  mTrig[2]=96201;
81  mTrig[3]=96211;
82  }
83  mDebug=debug;
84 
85  mBemcTables=new StBemcTables();
86  mEvent=new MyEvent();
87 }
88 
89 StPhotonMaker::~StPhotonMaker()
90 {
91  //destructor
92 }
93 
94 Int_t StPhotonMaker::Init()
95 {
96 
97  h_EvSum=new TH1F("h_EvSum","see code for details",22,-1.5,20.5);
98  h_bsmdeAdc=new TH1F("h_bsmdeAdc","adc",1200,-0.5,1199.5);
99  h_bsmdpAdc=new TH1F("h_bsmdpAdc","adc",1200,-0.5,1199.5);
100  h_btowAdc=new TH1F("h_btowAdc","adc",1200,-0.5,1199.5);
101  h_bsmdeEn=new TH1F("h_bsmdeEn","energy",55,-1.,10.);
102  h_btowEn=new TH1F("h_btowEn","energy",55,-1.,10.);
103  h_btowEnVsAdc=new TH2F("h_btowEnVsAdc","h_btowEnVsAdc",1200,-0.5,1199.5,100,0.,30.);
104 
105  mEventTree=new TTree("mEventTree","tree with my event structure");
106  mEventTree->Branch("branch","MyEvent",&mEvent);
107 
108  mN=0;
109  mRunId=0;
110  mRunPrev=0;
111  mEventId=0;
112  mTrigger=0;
113  mDate=0;
114  mTime=0;
115  mPs_mb=0;
116  mPs_mb2=0;
117  mPs_ht1=0;
118  mPs_ht2=0;
119 
120  return StMaker::Init();
121 }
122 
123 Int_t StPhotonMaker::Make()
124 {
125  mN++;
126  if(mDebug) cout<<endl<<"#### processing event "<<mN<<endl<<endl;
127  h_EvSum->Fill(1);
128 
129  StEmcGeom *bemcGeom=StEmcGeom::getEmcGeom("bemc");
130 
131  StEvent *event=0;
132  event = (StEvent*)GetInputDS("StEvent");
133  if(!event){
134  h_EvSum->Fill(2);
135  if(mDebug) cout << "************ StMyPhotonMaker::Make() no event pointer!! *********" << endl;
136  return kStOK;
137  }
138 
139  //test bbc coincidence:
140  //StL0Trigger *l0trig=(StL0Trigger*)event->l0Trigger();
141  //unsigned short dsmInput=l0trig->dsmInput();
142  //cout<<"testing dsm: "<<dsmInput<<endl;
143 
144 
145  if(!mHijing && !mMc && !mPythia && mAdcMaker->isCorrupted()){
146  if(mDebug) cout<<"event corrupted!"<<endl;
147  h_EvSum->Fill(-1);
148  return kStOK;
149  }
150 
151  mRunPrev=mRunId;
152  if(!mMc){
153  mRunId=event->runId();
154  mEventId=event->id();
155  }
156  else{
157  mEventId=mN;
158  mRunId=0;
159  }
160  if(mDebug) cout<<"starting run: "<<mRunId<<" and event: "<<mEventId<<endl;
161 
162 
163  if(mRunId!=mRunPrev&&mReal){
164  StDetectorDbTriggerID& DetDbTrigId = *(StDetectorDbTriggerID::instance());
165  for(UInt_t i=0;i<DetDbTrigId.getL0NumRows();i++){
166  if(mDebug) cout<<"prescales: "<<DetDbTrigId.getPsL0(i)<<endl<<endl;
167  if(DetDbTrigId.getL0OfflineTrgId(i)==mTrig[0]) mPs_mb=DetDbTrigId.getPsL0(i);
168  if(DetDbTrigId.getL0OfflineTrgId(i)==mTrig[1]) mPs_mb2=DetDbTrigId.getPsL0(i);
169  if(DetDbTrigId.getL0OfflineTrgId(i)==mTrig[2]) mPs_ht1=DetDbTrigId.getPsL0(i);
170  if(DetDbTrigId.getL0OfflineTrgId(i)==mTrig[3]) mPs_ht2=DetDbTrigId.getPsL0(i);
171  }
172  mBemcTables->loadTables((StMaker*)this);
173  }
174 
175  mDate=(Int_t)mDbMaker->GetDateTime().GetDate();
176  mTime=(Int_t)mDbMaker->GetDateTime().GetTime();
177 
178  //get trigger info
179  Float_t mips=-1.;
180  Float_t zdcEast=0.;
181  Float_t zdcWest=0.;
182  Float_t zdcVertex=0.;
183  Float_t bbcEast=0.;
184  Float_t bbcWest=0.;
185  StTriggerDetectorCollection *trigDetColl=event->triggerDetectorCollection();
186  if(trigDetColl){
187  StCtbTriggerDetector &CTB=trigDetColl->ctb();
188  mips=CTB.mips(0);
189  StZdcTriggerDetector &ZDC=trigDetColl->zdc();
190  zdcEast=ZDC.adcSum(east);
191  zdcWest=ZDC.adcSum(west);
192  zdcVertex=ZDC.vertexZ();
193  StBbcTriggerDetector &BBC=trigDetColl->bbc();
194  bbcEast=BBC.adcSumEast();
195  bbcWest=BBC.adcSumWest();
196  }
197 
198  //----- CREATE MyEvent -----
199  mEvent=new MyEvent(mRunId,mEventId,mDate,mTime);
200  if(!mEvent){
201  if(mDebug) cout<<"couldn't create MyEvent* !!"<<endl;
202  return kStErr;
203  }
204 
205  //extract MC data:
206  StPrimaryVertex *pVert=event->primaryVertex();
207  StMcVertex *mc_pVert=0;
208  StMcEvent* mc_event=0;
209  if(mPythia||mHijing){
210  mc_event=(StMcEvent*)GetInputDS("StMcEvent");
211  if(!mc_event){
212  h_EvSum->Fill(4);
213  if(mDebug) cout<<"StPhotonMaker::Make() no McEvent"<<endl;
214  return kStOK;
215  }
216  mc_pVert=mc_event->primaryVertex();
217  if(!mc_pVert){
218  h_EvSum->Fill(5);
219  if(mDebug) cout<<"StPhotonMaker::Make() no mc vertex!"<<endl;
220  return kStOK;
221  }
222  }
223  if(mMc||mEmbed){
224  mc_event=(StMcEvent*)GetInputDS("StMcEvent");
225  if(!mc_event){
226  h_EvSum->Fill(4);
227  if(mDebug) cout<<"StPhotonMaker::Make() no McEvent"<<endl;
228  return kStOK;
229  }
230  mc_pVert=mc_event->primaryVertex();
231  if(!mc_pVert){
232  h_EvSum->Fill(5);
233  if(mDebug) cout<<"StPhotonMaker::Make() no mc vertex!"<<endl;
234  return kStOK;
235  }
236  //get generated particles:
237  StPtrVecMcTrack& trcks=mc_pVert->daughters();
238  Int_t nDaughters=mc_pVert->numberOfDaughters();
239  if(mDebug) cout<<"number of MC part. from vertex: "<<nDaughters<<endl;
240  if(nDaughters==1){
241  StMcTrackIterator it=trcks.begin();
242  Int_t pid=(*it)->particleDefinition()->pdgEncoding();
243  if(mDebug) cout<<"daughter pid: "<<pid<<endl;
244  // pid=111 -> pi0; pid=221 -> eta; pid=22 -> gamma; pid=(-)2112 -> (anti-)proton; pid=130 -> kzero_L
245  // if we're talking mc gammas form the vertex, then we can fill one of the two branches
246  if(pid==22){
247  MyMcTrack *mytr=new MyMcTrack();
248  mytr->setId(pid);
249  mytr->setEnergy((*it)->energy());
250  if((*it)->stopVertex()){
251  if((*it)->stopVertex()->daughters().size()){
252  mytr->setStopRadius((*it)->stopVertex()->position().perp());
253  if(mDebug) cout<<"radius: "<<mytr->stopRadius()<<endl;
254  }
255  else if(mDebug) cout<<"no size"<<endl;
256  }
257  else if(mDebug) cout<<"no stop"<<endl;
258  mytr->setMomentum((*it)->momentum().x(),(*it)->momentum().y(),(*it)->momentum().z());
259 
260  StThreeVectorF v;
261  v.setX(bemcGeom->Radius()*cos((*it)->momentum().phi()));
262  v.setY(bemcGeom->Radius()*sin((*it)->momentum().phi()));
263  v.setZ(bemcGeom->Radius()/tan((*it)->momentum().theta()));
264  v=v + mc_pVert->position();
265 
266  mytr->setPosition(v.x(),v.y(),v.z());
267  //add generated mc particle to event
268  mEvent->setMcTrack(mytr);
269  delete mytr;
270  }
271  else if(pid==111||pid==221||pid==2112||pid==-2112||pid==211||pid==130){
272 
273  MyMcTrack *mytr=new MyMcTrack();
274  mytr->setId(pid);
275  mytr->setEnergy((*it)->energy());
276 
277  if((*it)->stopVertex()){
278  if((*it)->stopVertex()->daughters().size()){
279  mytr->setStopRadius((*it)->stopVertex()->position().perp());
280  }
281  }
282 
283  StThreeVectorF v;
284  v.setX(bemcGeom->Radius()*cos((*it)->momentum().phi()));
285  v.setY(bemcGeom->Radius()*sin((*it)->momentum().phi()));
286  v.setZ(bemcGeom->Radius()/tan((*it)->momentum().theta()));
287  v=v + mc_pVert->position();
288 
289  mytr->setMomentum((*it)->momentum().x(),(*it)->momentum().y(),(*it)->momentum().z());
290  mytr->setPosition(v.x(),v.y(),v.z());
291  //add generated mc particle to event
292  mEvent->setMcTrack(mytr);
293  delete mytr;
294 
295  //and the daughters -> photons!!
296  if((*it)->stopVertex()){
297  StMcVertex *sVert=(*it)->stopVertex();
298  StPtrVecMcTrack& gammaTracks=sVert->daughters();
299  for(StMcTrackIterator itd=gammaTracks.begin();itd!=gammaTracks.end();itd++){
300  Int_t pid2=(*itd)->particleDefinition()->pdgEncoding();
301  MyMcTrack *mytr2=new MyMcTrack();
302  mytr2->setId(pid2);
303  mytr2->setEnergy((*itd)->energy());
304 
305  if((*itd)->stopVertex()){
306  if((*itd)->stopVertex()->daughters().size()){
307  mytr2->setStopRadius((*it)->stopVertex()->position().perp());
308  }
309  }
310 
311  StThreeVectorF vv;
312  vv.setX(bemcGeom->Radius()*cos((*itd)->momentum().phi()));
313  vv.setY(bemcGeom->Radius()*sin((*itd)->momentum().phi()));
314  vv.setZ(bemcGeom->Radius()/tan((*itd)->momentum().theta()));
315  vv=vv + mc_pVert->position();
316 
317  mytr2->setMomentum((*itd)->momentum().x(),(*itd)->momentum().y(),(*itd)->momentum().z());
318  mytr2->setPosition(vv.x(),vv.y(),vv.z());
319 
320  mEvent->addMcPhoton(mytr2);
321  delete mytr2;
322  }
323  }
324  }
325  else
326  cout<<"what kind of MC particle is this??? (shouldn't see this)"<<endl;
327  }
328  else if(mMc && nDaughters>1){
329  if(mDebug) cout<<"more than one particle per event"<<endl;
330 
331  StMcTrackIterator it=trcks.begin();
332  Int_t pid=(*it)->particleDefinition()->pdgEncoding();
333  if(mDebug) cout<<"daughter pid: "<<pid<<endl;
334  if(pid==-2112 || pid==2112 || pid==221){
335  for(UInt_t i=0;i<trcks.size();i++){
336  Float_t neutronenergy=trcks[i]->energy();
337  MyMcTrack *myTR=new MyMcTrack();
338  myTR->setId(pid);
339  myTR->setEnergy(neutronenergy);
340  myTR->setMomentum(trcks[i]->momentum().x(),trcks[i]->momentum().y(),trcks[i]->momentum().z());
341  StThreeVectorF v;
342  v.setX(bemcGeom->Radius()*cos(trcks[i]->momentum().phi()));
343  v.setY(bemcGeom->Radius()*sin(trcks[i]->momentum().phi()));
344  v.setZ(bemcGeom->Radius()/tan(trcks[i]->momentum().theta()));
345  v=v + mc_pVert->position();
346  myTR->setPosition(v.x(),v.y(),v.z());
347  mEvent->addMcPion(myTR);
348  delete myTR;
349  }
350  }
351  }
352  }
353 
354  if(!pVert&&!mc_pVert){
355  h_EvSum->Fill(6);
356  if(mDebug) cout<<"StPhotonMaker::Make() no primary vertex"<<endl;
357  return kStOK;
358  }
359 
360  if(mPythia){
361  //get partonic pT (thanks Frank):
362  TDataSet *gEvent=GetDataSet("geant");
363  TDataSetIter geantDstI(gEvent);
364  //St_particle *particleTabPtr=(St_particle*)geantDstI("particle");
365  //particle_st* particleTable=particleTabPtr->GetTable();
366  //St_g2t_event *Pg2t_event=(St_g2t_event*)geantDstI("g2t_event");
367  //g2t_event_st *g2t_event1=Pg2t_event->GetTable();
368 
369  //Get parontic pT,xb1,xb2
370  St_g2t_pythia *Pg2t_pythia=(St_g2t_pythia*)geantDstI("g2t_pythia");
371  if(Pg2t_pythia){
372  g2t_pythia_st *g2t_pythia1=Pg2t_pythia->GetTable();
373  if(g2t_pythia1){
374  float hard_p=g2t_pythia1->hard_p;
375  mEvent->setPartonPt(hard_p);
376  }
377  }
378  }
379  if(mPythia||mHijing){
380  StPtrVecMcTrack& tracs=mc_event->tracks();
381  for(UInt_t i=0;i<tracs.size();i++){
382  //use geant id for pythia!
383  if(tracs[i]->geantId()==7){
384  Float_t pionenergy=tracs[i]->energy();
385  MyMcTrack *myTR=new MyMcTrack();
386  myTR->setId(111);
387  myTR->setEnergy(pionenergy);
388  myTR->setMomentum(tracs[i]->momentum().x(),tracs[i]->momentum().y(),tracs[i]->momentum().z());
389  StThreeVectorF v;
390  v.setX(bemcGeom->Radius()*cos(tracs[i]->momentum().phi()));
391  v.setY(bemcGeom->Radius()*sin(tracs[i]->momentum().phi()));
392  v.setZ(bemcGeom->Radius()/tan(tracs[i]->momentum().theta()));
393  v=v + mc_pVert->position();
394  myTR->setPosition(v.x(),v.y(),v.z());
395  mEvent->addMcPion(myTR);
396  delete myTR;
397  }
398  if(tracs[i]->geantId()==1){
399  Float_t photonenergy=tracs[i]->energy();
400  MyMcTrack *myTR=new MyMcTrack();
401  myTR->setId(22);
402 
403  //get parent id:
404  //int parentid=-1;
405  //if(tracs[i]->parent()){
406  // parentid=tracs[i]->parent()->geantId();
407  //}
408  //if(parentid==7) parentid=111;
409  //else if(parentid==17) parentid=221;
410  //else parentid=0;
411  //myTR->setParentId(parentid);
412  myTR->setEnergy(photonenergy);
413  myTR->setMomentum(tracs[i]->momentum().x(),tracs[i]->momentum().y(),tracs[i]->momentum().z());
414  StThreeVectorF v;
415  v.setX(bemcGeom->Radius()*cos(tracs[i]->momentum().phi()));
416  v.setY(bemcGeom->Radius()*sin(tracs[i]->momentum().phi()));
417  v.setZ(bemcGeom->Radius()/tan(tracs[i]->momentum().theta()));
418  v=v + mc_pVert->position();
419  myTR->setPosition(v.x(),v.y(),v.z());
420  mEvent->addMcPhoton(myTR);
421  delete myTR;
422  }
423  }
424  }
425 
426  if(!pVert&&mEmbed){
427  h_EvSum->Fill(7);
428  if(mDebug) cout<<"StPhotonMaker::Make() no primary vertex ****"<<endl;
429  return kStOK;
430  }
431  StThreeVectorF vPos,mc_vPos;
432  if(mMc||mEmbed||mPythia||mHijing){
433  mc_vPos=mc_pVert->position();
434  if(mMc||mPythia) vPos=mc_vPos;
435  }
436  if(mEmbed||mReal||mHijing){
437  if(pVert) vPos=pVert->position();
438  }
439 
440  //get trig per event
441  mTrigger=0;
442  if(mReal&&event->triggerIdCollection()&&event->triggerIdCollection()->nominal()){
443  if(event->triggerIdCollection()->nominal()->isTrigger(mTrig[0])||
444  event->triggerIdCollection()->nominal()->isTrigger(mTrig[1])) mTrigger+=1;
445  if(event->triggerIdCollection()->nominal()->isTrigger(mTrig[2])) mTrigger+=2;
446  if(event->triggerIdCollection()->nominal()->isTrigger(mTrig[3])) mTrigger+=4;
447  if(event->triggerIdCollection()->nominal()->isTrigger(2002)) mTrigger+=8;
448  }
449 
450  StEmcTriggerMaker *triggermaker=dynamic_cast<StEmcTriggerMaker*>(GetMaker("bemctrigger"));
451  assert(triggermaker);
452 
453  if(mMc||mEmbed){
454  mTrigger=1;
455  if(mDate<20060001){
456  if(triggermaker->is2005HT1()) mTrigger+=2;
457  if(triggermaker->is2005HT2()) mTrigger+=4;
458  }
459  else{
460  cout<<"check timestamp -> triggermaker"<<endl;
461  assert(0);
462  }
463  }
464  if(mPythia){
465  mTrigger=1;
466  //get result from triggerMaker, 2005 for now:
467  if(mDate<20050101) assert(0);
468  if(triggermaker->is2005HT1()) mTrigger+=2;
469  if(triggermaker->is2005HT2()) mTrigger+=4;
470  }
471  if(mHijing){
472  mTrigger=1;
473  }
474 
475  if(mTrigger==0){
476  if(mDebug) cout<<"nothing more useless than an event with no trigger!"<<endl;
477  return kStOK;
478  }
479 
480  StEmcCollection *emcCol=(StEmcCollection*)event->emcCollection();
481  if(!emcCol){
482  if(mDebug) cout<<"no emccollection!!"<<endl;
483  h_EvSum->Fill(8);
484  return kStOK;
485  }
486  StEmcDetector *emcDet=(StEmcDetector*)emcCol->detector(kBarrelEmcTowerId);
487  if(!emcDet){
488  if(mDebug) cout<<"no emcDet!!"<<endl;
489  h_EvSum->Fill(9);
490  return kStOK;
491  }
492  //get status of trigger tower
493  Int_t hiTowId=-1;
494  Int_t hiTowAdc=-1;
495  Int_t hiTowStatus=-1;
496  Float_t hiTowEnergy=0.;
497  if(mTrigger>1 || mEmbed || mMc){
498  for(UInt_t m=1;m<=emcDet->numberOfModules();m++){
499  StEmcModule *emcMod=emcDet->module(m);
500  if(!emcMod){
501  if(mDebug) cout<<"no emcMod for module "<<m<<endl;
502  continue;
503  }
504  StSPtrVecEmcRawHit& modHits=emcMod->hits();
505  for(UInt_t h=0;h<modHits.size();h++){
506  Int_t adc=modHits[h]->adc();
507  Float_t nrg=modHits[h]->energy();
508  //h_btowAdc->Fill(adc);
509  //h_btowEn->Fill(nrg);
510  //h_btowEnVsAdc->Fill(adc,nrg);
511  if(adc>hiTowAdc && nrg>0.){
512  hiTowAdc=adc;
513  hiTowEnergy=nrg;
514  bemcGeom->getId(m,modHits[h]->eta(),modHits[h]->sub(),hiTowId);
515  }
516  }
517  }
518  //get status
519  if(hiTowId>0) mBemcTables->getStatus(BTOW,hiTowId,hiTowStatus);
520  }
521 
522  StEmcDetector *bsmdeDet=(StEmcDetector*)emcCol->detector(kBarrelSmdEtaStripId);
523  if(!bsmdeDet){
524  return kStOK;
525  }
526  StEmcDetector *bsmdpDet=(StEmcDetector*)emcCol->detector(kBarrelSmdPhiStripId);
527  if(!bsmdpDet){
528  return kStOK;
529  }
530  for(UInt_t m2=1;m2<=bsmdeDet->numberOfModules();m2++){
531  StEmcModule *bsmdeMod=bsmdeDet->module(m2);
532  StSPtrVecEmcRawHit& bsmdeHits=bsmdeMod->hits();
533  for(UInt_t h2=0;h2<bsmdeHits.size();h2++){
534  h_bsmdeAdc->Fill(bsmdeHits[h2]->adc());
535  }
536  }
537  for(UInt_t m2=1;m2<=bsmdpDet->numberOfModules();m2++){
538  StEmcModule *bsmdpMod=bsmdpDet->module(m2);
539  StSPtrVecEmcRawHit& bsmdpHits=bsmdpMod->hits();
540  for(UInt_t h2=0;h2<bsmdpHits.size();h2++){
541  h_bsmdpAdc->Fill(bsmdpHits[h2]->adc());
542  }
543  }
544 
545 
546  StSPtrVecEmcPoint& barrelPoints=emcCol->barrelPoints();
547  Int_t nPrimTracks=0;
548  Int_t nGoodPrimaries=0;
549  Int_t nGoodPrimBarrel=0;
550  Int_t nGoodGlobals=0;
551  Float_t TpcPt=0.;
552  Float_t TpcPtBarrelWest=0.;
553  if(event->primaryVertex()){
554  nPrimTracks=pVert->numberOfDaughters();
555  StSPtrVecTrackNode& trackNode=event->trackNodes();
556  if(mDebug) cout<<"number of tracknodes: "<<trackNode.size()<<endl;
557  for(UInt_t j=0;j<trackNode.size();j++){
558  StPrimaryTrack* track=static_cast<StPrimaryTrack*>(trackNode[j]->track(primary));
559  if(track && track->flag()>0 && track->geometry()){
560  //get DCA
561  Double_t pathlength=track->geometry()->helix().pathLength(vPos,kFALSE);
562  StThreeVectorD v=track->geometry()->helix().at(pathlength) - vPos;
563  Float_t dca=v.mag();
564  //get number of fit points
565  const StTrackFitTraits& fitTraits=track->fitTraits();
566  unsigned short numFitPoints=fitTraits.numberOfFitPoints();
567  if(numFitPoints>15&&dca<3.){
568  nGoodPrimaries++;
569  Float_t tracEta=track->geometry()->momentum().pseudoRapidity();
570  if(tracEta>0.&&tracEta<1.){
571  nGoodPrimBarrel++;
572  TpcPtBarrelWest+=track->geometry()->momentum().perp();
573  }
574  TpcPt+=track->geometry()->momentum().perp();
575  }
576  }
577  StGlobalTrack* gtrack=static_cast<StGlobalTrack*>(trackNode[j]->track(global));
578  if(gtrack && gtrack->flag()>0 && gtrack->geometry()){
579  const StTrackFitTraits& fitTraitsG=gtrack->fitTraits();
580  unsigned short numFitPointsG=fitTraitsG.numberOfFitPoints();
581  if(numFitPointsG>15) nGoodGlobals++;
582  }
583  }
584  }
585 
586  //FTPC tracks for centrality, copied from "StRoot/StEventUtilities/StuFtpcRefMult.hh"
587  Int_t ftpcRefMultTracks=0;
588  if(event->primaryVertex()){
589  const StSPtrVecPrimaryTrack& tracks=pVert->daughters();
590  for(StSPtrVecPrimaryTrackConstIterator iter = tracks.begin(); iter != tracks.end(); iter++){
591  StTrack* track = (*iter);
592  if(!track->geometry()){
593  if(mDebug) cout<<"no track geometry! -> NEXT"<<endl;
594  continue;
595  }
596  //check if possible FTPC tracks
597  if(track->fitTraits().numberOfFitPoints()<6||(track->fitTraits().numberOfFitPoints()>11))
598  continue;
599  //check eta range and FTPC east
600  if(track->geometry()->momentum().pseudoRapidity()>-2.8||track->geometry()->momentum().pseudoRapidity()<=-3.8)
601  continue;
602  //check pt
603  if(track->geometry()->momentum().perp()>=3.)
604  continue;
605  //finally, check dca, if a track satisfies gets inside the if, count it.
606  StTrack *glt=track->node()->track(global);
607  if(!glt){
608  if(mDebug) cout<<"no global track!"<<endl;
609  continue;
610  }
611  if(!glt->geometry()){
612  if(mDebug) cout<<"no geom.!"<<endl;
613  continue;
614  }
615  if(glt->geometry()->helix().distance(vPos)<3.) ftpcRefMultTracks++;
616  }
617  }
618 
619  //start reco tree here:
620  StEmcPoint *pt=0;
621  Float_t EnergyNeutral=0.;
622  for(UInt_t i=0;i<barrelPoints.size();i++)
623  {
624  MyPoint *mypt=new MyPoint();
625  pt=(StEmcPoint*)barrelPoints[i];
626 
627  Float_t dist=9999.;
628  if(!calcDistanceTrackToPoint(pt,dist)){
629  cout<<"wrong bfield or no StEvent pointer (shouldn't see this)"<<endl;
630  }
631  if(pt->position().pseudoRapidity()>0.0)//only WEST
632  EnergyNeutral+=pt->energy();
633 
634  StPtrVecEmcCluster& etaClus=pt->cluster(kBarrelSmdEtaStripId);
635  StPtrVecEmcCluster& phiClus=pt->cluster(kBarrelSmdPhiStripId);
636  Int_t nEtaHits=0;
637  Int_t nPhiHits=0;
638  Float_t etaEnergy=0.;
639  Float_t phiEnergy=0.;
640  Float_t etaWidth=0.;
641  Float_t phiWidth=0.;
642  if(etaClus.size()>0){
643  StEmcCluster *clE=(StEmcCluster*)etaClus[0];
644  StPtrVecEmcRawHit& hE=clE->hit();
645  nEtaHits=hE.size();
646  etaEnergy=clE->energy();
647  etaWidth=clE->sigmaEta();
648  }
649  if(phiClus.size()>0){
650  StEmcCluster *clP=(StEmcCluster*)phiClus[0];
651  StPtrVecEmcRawHit& hP=clP->hit();
652  nPhiHits=hP.size();
653  phiEnergy=clP->energy();
654  phiWidth=clP->sigmaPhi();
655  }
656  mypt->setEnergy(pt->energy());
657  mypt->setQuality((Int_t)pt->quality());
658  mypt->setPosition(pt->position().x(),pt->position().y(),pt->position().z());
659  mypt->setDistanceToTrack(dist);
660  mypt->setHitsEta(nEtaHits);
661  mypt->setWidthEta(etaWidth);
662  mypt->setEnergyEta(etaEnergy);
663  mypt->setHitsPhi(nPhiHits);
664  mypt->setWidthPhi(phiWidth);
665  mypt->setEnergyPhi(phiEnergy);
666  mEvent->addPoint(mypt);
667 
668  delete mypt;
669  }
670 
671 
672  //fill MyEvent
673  mEvent->setPrescale(0,mPs_mb);
674  mEvent->setPrescale(1,mPs_mb2);
675  mEvent->setPrescale(2,mPs_ht1);
676  mEvent->setPrescale(3,mPs_ht2);
677  mEvent->setHighTowerAdc(hiTowAdc);
678  mEvent->setHighTowerId(hiTowId);
679  mEvent->setHighTowerStatus(hiTowStatus);
680  mEvent->setHighTowerEnergy(hiTowEnergy);
681  mEvent->setVertex(vPos.x(),vPos.y(),vPos.z());
682  mEvent->setTrigger(mTrigger);
683  mEvent->setGoodPrimaries(nGoodPrimaries);
684  mEvent->setGoodPrimBarrel(nGoodPrimBarrel);
685  mEvent->setGoodGlobals(nGoodGlobals);
686  mEvent->setRefMult(ftpcRefMultTracks);
687  mEvent->setEnergyInBarrel(EnergyNeutral);
688  mEvent->setMomentumInTpc(TpcPt);
689  mEvent->setMomentumInTpcWest(TpcPtBarrelWest);
690  mEvent->setCtbSum(mips);
691  mEvent->setBbcSumEast(bbcEast);
692  mEvent->setBbcSumWest(bbcWest);
693  mEvent->setZdcSumEast(zdcEast);
694  mEvent->setZdcSumWest(zdcWest);
695  mEvent->setZdcVertexZ(zdcVertex);
696  if(mHijing){
697  mEvent->setZdcVertexZ(mc_vPos.z());
698  }
699 
700  mEventTree->Fill();
701 
702  return kStOK;
703 }
704 
705 
706 
707 
708 
709 Int_t StPhotonMaker::Finish()
710 {
711  saveHistograms();
712  return kStOK;
713 }
714 
715 
716 void StPhotonMaker::saveHistograms()
717 {
718  if(mDebug) cout<<"************ saving histograms ****************"<<endl;
719  TFile *hfile=(TFile*)gROOT->FindObject(mFileName);
720  if(hfile) hfile->Close();
721  hfile=new TFile(mFileName,"RECREATE");
722  //hfile->SetCompressionLevel(1);
723  //----
724  h_EvSum->Write();
725  h_bsmdeAdc->Write();
726  h_bsmdpAdc->Write();
727  //h_bsmdeEn->Write();
728  //h_btowAdc->Write();
729  //h_btowEn->Write();
730  //h_btowEnVsAdc->Write();
731 
732  mEventTree->Write();
733 
734  //----
735  hfile->Close();
736 }
737 
738 void StPhotonMaker::setDbMaker(St_db_Maker *dbM)
739 {
740  mDbMaker=dbM;
741 }
742 void StPhotonMaker::setAdcMaker(StEmcADCtoEMaker *adc)
743 {
744  mAdcMaker=adc;
745 }
746 Bool_t StPhotonMaker::calcDistanceTrackToPoint(StEmcPoint *point,Float_t &distanceToTrack)
747 {
748  StEvent *event = (StEvent*)this->GetInputDS("StEvent");
749  if(!event){
750  if(mDebug) cout<<"can't get Event pointer"<<endl;
751  return kFALSE;
752  }
753 
754  Double_t bFld=0.;
755  StEventSummary* summary = event->summary();
756  if(summary){
757  bFld = summary->magneticField()/10.; // bFld in Tesla
758  if(mDebug) cout<<"summary()->magneticField()="<<bFld<<"(Tesla)"<<endl;
759  }
760  if(TMath::Abs(bFld)<0.01&&!mMc){
761  if(mDebug) cout<<"wrong mBField!"<<endl;
762  return kFALSE;
763  }
764  StEmcPosition* emcPosition = new StEmcPosition();
765  StThreeVectorD pos, mom;
766  StSPtrVecTrackNode& trackNodes = event->trackNodes();
767  StTrack* track;
768  distanceToTrack=999.;
769  for (size_t nodeIndex=0; nodeIndex<trackNodes.size();nodeIndex++){
770  size_t numberOfTracksInNode=trackNodes[nodeIndex]->entries(global);
771  for(size_t trackIndex=0; trackIndex<numberOfTracksInNode;trackIndex++){
772  track = trackNodes[nodeIndex]->track(primary,trackIndex);
773  if (track && track->flag()>=0){
774  if(track->geometry()->momentum().mag()<.5) continue;
775  StPhysicalHelixD helix=track->outerGeometry()->helix();
776  Bool_t projOK=emcPosition->projTrack(&pos,&mom,track,bFld,230.705,1);
777  if(!projOK) continue;
778  Float_t d=(pos - point->position()).mag();//in cm
779  if(d<distanceToTrack&&d>0.) distanceToTrack=d;
780  }
781  }
782  }
783  delete emcPosition;
784  return kTRUE;
785 }
St_db_Maker
Definition: St_db_Maker.h:141
StPhotonMaker::Finish
virtual Int_t Finish()
Definition: StPhotonMaker.cxx:709
StEmcADCtoEMaker::isCorrupted
Bool_t isCorrupted()
Returns if BTOW is corrupted or not.
Definition: StEmcADCtoEMaker.cxx:239
StHelix::distance
double distance(const StThreeVector< double > &p, bool scanPeriods=true) const
minimal distance between point and helix
Definition: StHelix.cc:240
MyPoint
Definition: MyPoint.h:7
StHelix::pathLength
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
Definition: StHelix.cc:351
StDetectorDbTriggerID
Accessor to the database for trigger id information.
Definition: StDetectorDbTriggerID.h:47
StBemcTables::loadTables
void loadTables(StMaker *anyMaker)
load tables.
Definition: StBemcTables.cxx:73
kStOK
Definition: Stypes.h:40
StPhotonMaker::Make
virtual Int_t Make()
Definition: StPhotonMaker.cxx:123
StBemcTables::getStatus
void getStatus(Int_t det, Int_t softId, Int_t &status, const char *option="") const
Return status.
Definition: StBemcTables.cxx:341
StEmcPosition::projTrack
Bool_t projTrack(StThreeVectorD *atFinal, StThreeVectorD *momentumAtFinal, const StTrack *const track, Double_t magField, Double_t radius=225.405, Int_t option=1) const
Track projection utility.
Definition: StEmcPosition.cxx:139
kStErr
Definition: Stypes.h:44
StMcEvent
Event data structure to hold all information from a Monte Carlo simulation. This class is the interfa...
Definition: StMcEvent.hh:169
Generated by   doxygen 1.8.5