StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
WeventDisplay.cxx
1 // $Id: WeventDisplay.cxx,v 1.6 2010/01/10 03:01:37 balewski Exp $
2 //
3 //*-- Author : Jan Balewski, MIT
4 
5 
6 #include <TH1.h>
7 #include <TH2.h>
8 #include <TMath.h>
9 #include <TCanvas.h>
10 #include <TStyle.h>
11 #include <TFile.h>
12 #include <TText.h>
13 #include <TLine.h>
14 #include <TList.h>
15 #include <TBox.h>
16 #include <TPaveText.h>
17 #include <TPad.h>
18 #include <TEllipse.h>
19 #include <stdio.h>
20 
21 //MuDst
22 #include <StMuDSTMaker/COMMON/StMuDstMaker.h>
23 #include <StMuDSTMaker/COMMON/StMuDst.h>
24 #include <StMuDSTMaker/COMMON/StMuTriggerIdCollection.h>
25 #include <StMuDSTMaker/COMMON/StMuEvent.h>
26 #include <StMuDSTMaker/COMMON/StMuTrack.h>
27 #include <StMuDSTMaker/COMMON/StMuPrimaryVertex.h>
28 #include "StEmcUtil/geometry/StEmcGeom.h"
29 
30 #include "St2009WMaker.h"
31 #include "WanaConst.h"
32 #include "WeventDisplay.h"
33 //-----------------------------
34 //-----------------------------
35 WeventDisplay::WeventDisplay( St2009WMaker* mk, int mxEv) {
36  maxEve=mxEv;
37  wMK=mk;
38  const float PI=TMath::Pi();
39  const char cPlane[ mxBSmd]={'E','P'};
40  char txt1[100], txt2[1000];
41 
42  // barrel
43  etaBL_ln=new TLine(-1,-3.2,1,3.2); etaBL_ln->SetLineColor(kBlue);
44  etaBR_ln=new TLine(-1,-3.2,1,3.2); etaBR_ln->SetLineColor(kBlue);
45 
46  etaEL_ln=new TLine(2,-3.2,3,3.2); etaEL_ln->SetLineColor(kGreen);
47  etaER_ln=new TLine(2,-3.2,3,3.2); etaER_ln->SetLineColor(kGreen);
48 
49  bxT=new TBox(-1.3,0, 1.3,1.); bxT->SetFillStyle(0); bxT->SetLineStyle(2);
50  bxE=new TBox(-1.3,0, 1,2.); bxE->SetFillStyle(0); bxE->SetLineStyle(2);
51  if(wMK->par_useEtow>=2) bxE->SetX2(2.2);
52 
53 
54  hEmcET=new TH2F("eveBtowET","EMC ET sum, Z=[0.3,30]GeV; event eta ; phi",38,-1.4,2.4,63,-PI,PI);
55 
56  hTpcET=new TH2F("eveTpcET","TPC PT sum, Z[0.3,10]GeV/c; event eta ; phi",28,-1.4,1.4,63,-PI,PI);
57 
58  for(int iep=0;iep<mxBSmd;iep++) {
59  sprintf(txt1,"eveBsmdAdc_%c",cPlane[iep]);
60  sprintf(txt2,"BSMD_%c ADC sum Z=[30,1000]; event eta ; phi",cPlane[iep]);
61  hBsmdAdc[iep]=new TH2F(txt1,txt2,26,-1.3,1.3,63,-PI,PI);
62  }
63 
64 
65 #if 0
66  //---- smart code from Willie
67  float etabinsA[1+mxBetaStrMod*2], etaphibinsA[mxBMod2Pi+1];
68  for(int i=0;i<mxBetaStrMod+1;i++)
69  etabinsA[mxBetaStrMod-i]=-(etabinsA[mxBetaStrMod+i]=wMK->mSmdEGeom->EtaB()[i]);
70  for(int i=0;i<mxBMod2Pi+1;i++)
71  etaphibinsA[i]=wMK->mSmdEGeom->PhiB()[i];
72  hBsmdEtaAdc=new TH2F("eveBsmdEtaAdc"," Event: BSMD-Eta ADC vs. eta & phi; pseudorapidity; azimuth",mxBetaStrMod*2,etabinsA,mxBMod2Pi,etaphibinsA);
73 #endif
74 
75 }
76 
77 
78 
79 //-----------------------------
80 //-----------------------------
81 void
82 WeventDisplay::clear(){
83  hEmcET->Reset();
84  hTpcET->Reset();
85  for(int iep=0;iep<mxBSmd;iep++) hBsmdAdc[iep]->Reset();
86 
87 }
88 
89 //-----------------------------
90 //-----------------------------
91 void
92 WeventDisplay::draw( const char *tit,int eveID, int daqSeq, int runNo, WeveVertex myV, WeveEleTrack myTr){
93  if(maxEve<=0) return;
94  maxEve--;
95  TStyle *myStyle=new TStyle();
96  myStyle->cd();
97  myStyle->SetPalette(1,0);
98  myStyle->SetOptStat(1000010);
99 
100  char txt[1000];
101  sprintf(txt,"display-%s_run%d.eventId%06dvert%d",tit,runNo,eveID,myV.id);
102 
103  printf("WeventDisplay::Draw %s\n",txt);
104  TCanvas *c0=new TCanvas(txt,txt,850,600);
105  c0->cd();
106 
107  TString tt=txt;
108  TPad *cU = new TPad(tt+"U", tt+"U",0.,0.2,1.,1.); cU->Draw();
109  TPad *cD = new TPad(tt+"D", tt+"D",0.,0.,1.,0.2); cD->Draw();
110 
111  cU->cd();
112  TPad *cU1 = new TPad(tt+"U1", tt+"U1",0.,0.,0.24,1.); cU1->Draw();
113  TPad *cU2 = new TPad(tt+"U2", tt+"U2",0.24,0.,0.55,1.); cU2->Draw();
114  TPad *cU3 = new TPad(tt+"U3", tt+"U3",0.55,0.,1.,1.); cU3->Draw();
115  cU3->Divide(2,1);
116 
117 
118  cU1->cd(); hTpcET->Draw("colz");
119 
120  TVector3 rW=myTr.pointTower.R;
121  rW.SetZ(rW.z()-myV.z);
122 
123  TEllipse *te1=new TEllipse(rW.Eta(),rW.Phi(), 0.1,0.1);
124  te1->SetFillStyle(0);te1->SetLineStyle(3); te1->SetLineColor(kMagenta);
125  TEllipse *te2=new TEllipse(rW.Eta(),rW.Phi(), 0.7, 0.7);
126  te2->SetFillStyle(0);te2->SetLineStyle(3); te2->SetLineColor(kBlack);
127 
128  TVector3 rA=-rW; // away direction
129  bxT->SetY1(rA.Phi() - wMK->par_awayDeltaPhi);
130  bxT->SetY2(rA.Phi() + wMK->par_awayDeltaPhi);
131 
132  bxE->SetY1(rA.Phi() - wMK->par_awayDeltaPhi);
133  bxE->SetY2(rA.Phi() + wMK->par_awayDeltaPhi);
134 
135 
136  te1->Draw(); te2->Draw(); bxT->Draw("l");
137  etaBL_ln->Draw(); etaBR_ln->Draw(); etaEL_ln->Draw();
138 
139  cU2->cd(); hEmcET->Draw("colz");
140  te1->Draw(); te2->Draw(); bxE->Draw("l");
141  etaBL_ln->Draw(); etaBR_ln->Draw();
142  etaEL_ln->Draw(); etaER_ln->Draw();
143 
144  for(int iep=0;iep<mxBSmd;iep++) {
145  cU3->cd(1+iep);
146  hBsmdAdc[iep]->Draw("colz");
147  te1->Draw(); te2->Draw(); bxT->Draw("l");
148  etaBL_ln->Draw(); etaBR_ln->Draw();
149  }
150 
151  //........... text information .............
152  TPaveText *pvt = new TPaveText(0,0.,1,1,"br");
153  cD->cd();
154  sprintf(txt,"run=%d eveID=%05d daq=%d vertex:ID=%d Z=%.0fcm",runNo,eveID,daqSeq,myV.id,myV.z);
155  printf("WeventDisplay::Event ID %s\n",txt);
156  pvt->AddText(txt);
157 
158  sprintf(txt,"TPC PT(GeV/c) near=%.1f away=%.1f ", myTr.nearTpcPT, myTr.awayTpcPT);
159  printf("WeventDisplay::Event TPC %s\n",txt);
160  pvt->AddText(txt);
161 
162  sprintf(txt,"BTOW ET/GeV: 2x2=%.1f near= %.1f away= %.1f",myTr.cluster.ET,myTr.nearBtowET,myTr.awayBtowET);
163  printf("WeventDisplay:: BTOW %s\n",txt);
164  pvt->AddText(txt);
165 
166  sprintf(txt,"Emc (Btow+Etow) ET/GeV: near= %.1f away= %.1f",myTr.nearEmcET,myTr.awayEmcET);
167  printf("WeventDisplay:: BTOW+ETOW %s\n",txt);
168  pvt->AddText(txt);
169 
170  sprintf(txt,"total ET/GeV: near= %.1f away= %.1f ptBalance= %.1f",myTr.nearTotET,myTr.awayTotET,myTr.ptBalance.Perp());
171  printf("WeventDisplay:: BTOW %s\n",txt);
172  pvt->AddText(txt);
173 
174 
175  pvt->Draw();
176 
177  c0->Print();
178  // save dump of histos
179  sprintf(txt,"display-%s_run%d.eventId%05dvert%d.root",tit,runNo,eveID,myV.id);
180  TFile hf(txt,"recreate");
181  if(hf.IsOpen()) {
182  hEmcET->Write();
183  hTpcET->Write();
184  for(int iep=0;iep<mxBSmd;iep++) hBsmdAdc[iep]->Write();
185  hf.Close();
186  }
187 }
188 
189 //-----------------------------
190 //-----------------------------
191 void
192 WeventDisplay::exportEvent( const char *tit, WeveVertex myV, WeveEleTrack myTr){
193  if(maxEve<=0) return;
194  clear();
195  int eveId=wMK->mMuDstMaker->muDst()->event()->eventId();
196  int runNo=wMK->mMuDstMaker->muDst()->event()->runId();
197  const char *afile = wMK->mMuDstMaker->GetFile();
198  int len=strlen(afile);
199  int daqSeq=atoi(afile+(len-18));
200  // printf("DDD %s len=%d %d =%s=\n",afile,len,daqSeq,afile+(len-15));
201 
202  TVector3 rTw=myTr.cluster.position;
203  rTw.SetZ(rTw.z()-myV.z);
204  printf("#xcheck-%s run=%d daqSeq=%d eveID=%7d vertID=%2d zVert=%.1f prTrID=%4d prTrEta=%.3f prTrPhi/deg=%.1f globPT=%.1f hitTwId=%4d twAdc=%.1f clEta=%.3f clPhi/deg=%.1f clET=%.1f\n",tit,
205  runNo,daqSeq,eveId,myV.id,myV.z,
206  myTr.prMuTrack->id(),myTr.prMuTrack->eta(),myTr.prMuTrack->phi()/3.1416*180.,myTr.glMuTrack->pt(),
207  myTr.pointTower.id,wMK->wEve.bemc.adcTile[kBTow][myTr.pointTower.id-1],
208  rTw.Eta(),rTw.Phi()/3.1416*180.,myTr.cluster.ET);
209 
210  float zVert=myV.z;
211  printf("WeventDisplay-%s::export run=%d eve=%d\n",tit,runNo,eveId);
212 
213  //.... process BTOW hits
214  for(int i=0;i< mxBtow;i++) {
215  float ene=wMK->wEve.bemc.eneTile[kBTow][i];
216  if(ene<=0) continue;
217  TVector3 primP=wMK->positionBtow[i]-TVector3(0,0,zVert);
218  primP.SetMag(ene); // it is 3D momentum in the event ref frame
219  float ET=primP.Perp();
220 
221  float eveEta=primP.Eta();
222  float evePhi=primP.Phi();
223  hEmcET->Fill(eveEta,evePhi,ET);
224  }
225 
226  int par_useEtow = wMK->par_useEtow; //get ETOW flag
227  //.... store ETOW hits
228  if(par_useEtow >= 1) {
229  for(int i=0;i<mxEtowPhiBin;i++){
230  for(int j=0;j<mxEtowEta;j++){
231  float ene=wMK->wEve.etow.ene[i][j];
232  if(ene<=0) continue;
233  TVector3 primP=wMK->positionEtow[i][j]-TVector3(0,0,zVert);
234  primP.SetMag(ene); // it is 3D momentum in the event ref frame
235  float ET=primP.Perp();
236 
237  float eveEta=primP.Eta();
238  float evePhi=primP.Phi();
239  hEmcET->Fill(eveEta,evePhi,ET);
240  }
241  }
242  }
243 
244  hEmcET->SetMinimum(0.3); hEmcET->SetMaximum(30.);
245  // compute approximate event eta for barrel
246  float x,y,z;
247  float Rcylinder= wMK->mBtowGeom->Radius();
248  assert(wMK->mBtowGeom->getXYZ(20,x,y,z)==0); // this is approximate Z of last tower
249  TVector3 rL(Rcylinder,0,z+myV.z);
250  TVector3 rR(Rcylinder,0,z-myV.z);
251  float etaL=-rL.Eta(), etaR=rR.Eta();
252  etaBL_ln->SetX1(etaL); etaBL_ln->SetX2(etaL);
253  etaBR_ln->SetX1(etaR); etaBR_ln->SetX2(etaR);
254  if(wMK->par_useEtow<2) bxE->SetX2(etaR+0.05); // skip endcap
255 
256  // compute approximate event eta for Endcap
257  rL=TVector3(0,214,270-myV.z); // detector eta~1.06
258  rR=TVector3(0,77,270-myV.z); // detector eta 2.0
259  etaL=rL.Eta(); etaR=rR.Eta();
260  etaEL_ln->SetX1(etaL); etaEL_ln->SetX2(etaL);
261  etaER_ln->SetX1(etaR); etaER_ln->SetX2(etaR);
262 
263 
264  //... TPC
265  hTpcET->SetMinimum(0.3);hTpcET->SetMaximum(10.);
266  getPrimTracks( myV.id);
267 
268  //.... BSMD-Eta, -Phi
269 
270  for(int iep=0;iep<mxBSmd;iep++) {
271  hBsmdAdc[iep]->SetMinimum(30);hBsmdAdc[iep]->SetMaximum(999);
272  for(int i=0;i< mxBStrips;i++) {
273  float adc=wMK->wEve.bemc.adcBsmd[iep][i];
274  if(adc<=0) continue;
275  TVector3 r=wMK->positionBsmd[iep][i];
276  float z1=r.z()-zVert;
277  r.SetZ(z1);
278  hBsmdAdc[iep]->Fill(r.Eta(),r.Phi(),adc);
279  }
280  }// end of eta,phi-planes
281 
282  //.... produce plot & save
283  draw(tit,eveId, daqSeq,runNo,myV, myTr);
284  export2sketchup(tit,myV, myTr);
285 }
286 
287 //-----------------------------
288 //-----------------------------
289 void
290 WeventDisplay::getPrimTracks( int vertID) {
291  assert(vertID>=0);
292  assert(vertID<(int)wMK->mMuDstMaker->muDst()->numberOfPrimaryVertices());
293  StMuPrimaryVertex* V=wMK-> mMuDstMaker->muDst()->primaryVertex(vertID);
294  assert(V);
295  wMK-> mMuDstMaker->muDst()->setVertexIndex(vertID);
296  float rank=V->ranking();
297  assert(rank>0);
298 
299  Int_t nPrimTrAll=wMK->mMuDstMaker->muDst()->GetNPrimaryTrack();
300  for(int itr=0;itr<nPrimTrAll;itr++) {
301  StMuTrack *prTr=wMK->mMuDstMaker->muDst()->primaryTracks(itr);
302  if(prTr->flag()<=0) continue;
303  if(prTr->flag()!=301) continue;// TPC-only regular tracks
304  float hitFrac=1.*prTr->nHitsFit()/prTr->nHitsPoss();
305  if(hitFrac<wMK->par_nHitFrac) continue;
306  StThreeVectorF prPvect=prTr->p();
307  TVector3 primP=TVector3(prPvect.x(),prPvect.y(),prPvect.z());
308  float pT=prTr->pt();
309  hTpcET->Fill(prTr->eta(),prTr->phi(),pT);
310 
311  }
312 }
313 
314 
315 //-----------------------------
316 //-----------------------------
317 void
318 WeventDisplay::export2sketchup( const char *tit, WeveVertex myV, WeveEleTrack myTr){
319  int eveId=wMK->mMuDstMaker->muDst()->event()->eventId();
320  int runNo=wMK->mMuDstMaker->muDst()->event()->runId();
321  char txt[1000];
322  sprintf(txt,"display3D-%s_run%d.eventId%05d_vert%d.txt",tit,runNo,eveId,myV.id);
323  FILE *fd=fopen(txt,"w"); assert(fd);
324 
325  //........ DUMP PRIM TRACKS..........
326  int vertID=myV.id;
327  assert(vertID>=0);
328  assert(vertID<(int)wMK->mMuDstMaker->muDst()->numberOfPrimaryVertices());
329  StMuPrimaryVertex* V=wMK-> mMuDstMaker->muDst()->primaryVertex(vertID);
330  assert(V);
331  wMK-> mMuDstMaker->muDst()->setVertexIndex(vertID);
332  float rank=V->ranking();
333  assert(rank>0);
334  const StThreeVectorF &rV=V->position();
335  Int_t nPrimTrAll=wMK->mMuDstMaker->muDst()->GetNPrimaryTrack();
336  for(int itr=0;itr<nPrimTrAll;itr++) {
337  StMuTrack *prTr=wMK->mMuDstMaker->muDst()->primaryTracks(itr);
338  if(prTr->flag()<=0) continue;
339  if(prTr->flag()!=301) continue;// TPC-only regular tracks
340  float hitFrac=1.*prTr->nHitsFit()/prTr->nHitsPoss();
341  if(hitFrac<wMK->par_nHitFrac) continue;
342  prTr->p();
343  fprintf(fd,"track V %.1f %.3f %.3f primP:PT:eta:phi:Q %.1f %.3f %.3f %d\n",rV.x(),rV.y(),rV.z(), prTr->p().perp(),prTr->p().pseudoRapidity(),prTr->p().phi(),prTr->charge());
344  }
345 
346  //........DUMP BTOW towers
347  float Rcylinder= wMK->mBtowGeom->Radius(), Rcylinder2=Rcylinder*Rcylinder;
348  for(int i=0;i< mxBtow;i++) {
349  float ene=wMK->wEve.bemc.eneTile[kBTow][i];
350  if(ene<=0) continue;
351  float delZ=wMK->positionBtow[i].z()-myV.z;
352  float e2et=Rcylinder/sqrt(Rcylinder2+delZ*delZ);
353  float ET=ene*e2et;
354  float detEta=wMK->positionBtow[i].Eta();
355  float detPhi=wMK->positionBtow[i].Phi();
356  fprintf(fd,"btow V %.1f %.3f %.3f eveET:detEta:detPhi %.3f %.3f %.3f\n",rV.x(),rV.y(),rV.z(),ET, detEta,detPhi);
357  }
358 
359  const char cPlane[ mxBSmd]={'E','P'};
360  //........DUMP BSMD hits
361  for(int iep=0;iep<mxBSmd;iep++) {
362  for(int i=0;i< mxBStrips;i++) {
363  float adc=wMK->wEve.bemc.adcBsmd[iep][i];
364  if(adc<=0) continue;
365  TVector3 r=wMK->positionBsmd[iep][i];
366  fprintf(fd,"bsmd%c V %.1f %.3f %.3f adc:detEta:detPhi %.3f %.3f %.3f\n",cPlane[iep],rV.x(),rV.y(),rV.z(),adc,r.Eta(),r.Phi() );
367  }
368  }
369 
370  //........DUMP ETOW towers
371  for(int iphi=0; iphi<mxEtowPhiBin; iphi++){
372  for(int ieta=0; ieta<mxEtowEta; ieta++){//sum all eta rings
373  float ene=wMK->wEve.etow.ene[iphi][ieta];
374  if(ene<=0) continue; //skip towers with no energy
375  TVector3 detP=wMK->positionEtow[iphi][ieta];
376  TVector3 primP=detP-TVector3(0,0,myV.z);
377  primP.SetMag(ene); // it is 3D momentum in the event ref frame
378  fprintf(fd,"etow V %.1f %.3f %.3f eveET:detEta:detPhi %.3f %.3f %.3f\n",rV.x(),rV.y(),rV.z(),primP.Perp(), detP.Eta(),detP.Phi());
379  }
380  }
381 
382  //.... DUMP reco electron ...
383  float eleET=myTr.cluster.ET;
384  fprintf(fd,"recoElectron V %.1f %.3f %.3f ET:detEta:detPhi %.3f %.3f %.3f\n",rV.x(),rV.y(),rV.z() , eleET,myTr.pointTower.R.Eta(),myTr.pointTower.R.Phi());
385 
386 
387  //... close event file
388  fclose(fd);
389 
390 
391 }
392 
393 
394 // $Log: WeventDisplay.cxx,v $
395 // Revision 1.6 2010/01/10 03:01:37 balewski
396 // cleanup & nicer histos
397 //
398 // Revision 1.5 2010/01/09 00:07:16 stevens4
399 // add jet finder
400 //
401 // Revision 1.4 2010/01/06 19:16:48 stevens4
402 // track cuts now on primary component, cleanup
403 //
404 // Revision 1.3 2010/01/06 04:22:15 balewski
405 // added Q/PT plot for Zs, more cleanup
406 //
407 // Revision 1.2 2009/12/10 16:01:31 stevens4
408 // fixed zero-length vector
409 //
410 // Revision 1.1 2009/11/23 23:00:18 balewski
411 // code moved spin-pool
412 //
StMuDst::primaryVertex
static StMuPrimaryVertex * primaryVertex()
return pointer to current primary vertex
Definition: StMuDst.h:322
StMuDstMaker::muDst
StMuDst * muDst()
Definition: StMuDstMaker.h:425
StMuTrack::pt
Double_t pt() const
Returns pT at point of dca to primary vertex.
Definition: StMuTrack.h:256
StMuTrack::id
short id() const
Returns the track id(or key), is unique for a track node, i.e. global and primary tracks have the sam...
Definition: StMuTrack.h:228
St2009WMaker
muDst based extraction of W-signal from pp500 data from 2009
Definition: St2009WMaker.h:50
StMuDstMaker::GetFile
virtual const char * GetFile() const
Returns name of current input or output file, depending on mode (GetFileName does the same...
Definition: StMuDstMaker.cxx:762
StMuTrack::nHitsFit
UShort_t nHitsFit() const
Return total number of hits used in fit.
Definition: StMuTrack.h:239
StMuTrack::flag
short flag() const
Returns flag, (see StEvent manual for type information)
Definition: StMuTrack.h:230
StMuTrack::charge
Short_t charge() const
Returns charge.
Definition: StMuTrack.h:255
StMuTrack::p
const StThreeVectorF & p() const
Returns 3-momentum at dca to primary vertex.
Definition: StMuTrack.h:259
StMuTrack::eta
Double_t eta() const
Returns pseudo rapidity at point of dca to primary vertex.
Definition: StMuTrack.h:257
StMuDst::primaryTracks
static TObjArray * primaryTracks()
returns pointer to a list of tracks belonging to the selected primary vertex
Definition: StMuDst.h:301
StMuTrack::nHitsPoss
UShort_t nHitsPoss() const
Return number of possible hits on track.
Definition: StMuTrack.cxx:242
StMuDst::event
static StMuEvent * event()
returns pointer to current StMuEvent (class holding the event wise information, e.g. event number, run number)
Definition: StMuDst.h:320
StMuTrack::phi
Double_t phi() const
Returns phi at point of dca to primary vertex.
Definition: StMuTrack.h:258
StMuDst::setVertexIndex
static void setVertexIndex(Int_t vtx_id)
Set the index number of the current primary vertex (used by both primaryTracks() functions and for St...
Definition: StMuDst.cxx:273
Generated by   doxygen 1.8.5