StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
WeventDisplay.cxx
1 // $Id: WeventDisplay.cxx,v 1.6 2012/09/21 16:59:10 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 <TLatex.h>
14 #include <TLine.h>
15 #include <TList.h>
16 #include <TBox.h>
17 #include <TPaveText.h>
18 #include <TPad.h>
19 #include <TEllipse.h>
20 #include <stdio.h>
21 
22 //MuDst
23 #include <StMuDSTMaker/COMMON/StMuDstMaker.h>
24 #include <StMuDSTMaker/COMMON/StMuDst.h>
25 #include <StMuDSTMaker/COMMON/StMuTriggerIdCollection.h>
26 #include <StMuDSTMaker/COMMON/StMuEvent.h>
27 #include <StMuDSTMaker/COMMON/StMuTrack.h>
28 #include <StMuDSTMaker/COMMON/StMuPrimaryVertex.h>
29 #include "StEmcUtil/geometry/StEmcGeom.h"
30 
31 //Esmd geometry for plots
32 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
33 #include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h"
34 
35 #include "St2011WMaker.h"
36 #include "WanaConst.h"
37 #include "WeventDisplay.h"
38 //-----------------------------
39 //-----------------------------
40 WeventDisplay::WeventDisplay( St2011WMaker* mk, int mxEv) {
41  maxEve=mxEv;
42  wMK=mk;
43  const float PI=TMath::Pi();
44  const char cPlane[ mxBSmd]={'E','P'};
45  const char cEsmdPlane[ mxEsmdPlane]={'U','V'};
46  char txt1[100], txt2[1000];
47 
48  // barrel
49  etaBL_ln=new TLine(-1,-3.2,1,3.2); etaBL_ln->SetLineColor(kBlue);
50  etaBR_ln=new TLine(-1,-3.2,1,3.2); etaBR_ln->SetLineColor(kBlue);
51 
52  etaEL_ln=new TLine(2,-3.2,3,3.2); etaEL_ln->SetLineColor(kGreen);
53  etaER_ln=new TLine(2,-3.2,3,3.2); etaER_ln->SetLineColor(kGreen);
54 
55  bxT=new TBox(-1.3,0, 1.3,1.); bxT->SetFillStyle(0); bxT->SetLineStyle(2);
56  bxE=new TBox(-1.3,0, 1,2.); bxE->SetFillStyle(0); bxE->SetLineStyle(2);
57  bxE->SetX2(2.2);
58 
59 
60  hEmcET=new TH2F("eveBtowET","EMC ET sum, Z=[0.3,30]GeV; event eta ; phi",38,-1.4,2.4,63,-PI,PI);
61 
62  hTpcET=new TH2F("eveTpcET","TPC PT sum, Z[0.3,10]GeV/c; event eta ; phi",32,-1.6,1.6,63,-PI,PI);
63 
64  for(int iep=0;iep<mxBSmd;iep++) {
65  sprintf(txt1,"eveBsmdAdc_%c",cPlane[iep]);
66  sprintf(txt2,"BSMD_%c ADC sum Z=[30,1000]; event eta ; phi",cPlane[iep]);
67  hBsmdAdc[iep]=new TH2F(txt1,txt2,26,-1.3,1.3,63,-PI,PI);
68  }
69 
70  //ESMD shower shape
71  for(int iuv=0; iuv<mxEsmdPlane; iuv++){
72  sprintf(txt1,"eveEsmdShower_%c",cEsmdPlane[iuv]);
73  sprintf(txt2,"ESMD_%c Shower Shape; i_strip position - track position (cm) ; MeV",cEsmdPlane[iuv]);
74  hEsmdShower[iuv]=new TH1F(txt1,txt2,41,-10.25,10.25);
75  bxEs[iuv]=new TBox(-5,0, 3,20.); bxEs[iuv]->SetFillStyle(3002); bxEs[iuv]->SetFillColor(kYellow);
76 
77  }
78  hEsmdXpt=new TH2F("eveEsmdXpt","ESMD Cross Point; X (cm); Y (cm)",100,0.,100,100,0.,100.);
79 
80 
81 #if 0
82  //---- smart code from Willie
83  float etabinsA[1+mxBetaStrMod*2], etaphibinsA[mxBMod2Pi+1];
84  for(int i=0;i<mxBetaStrMod+1;i++)
85  etabinsA[mxBetaStrMod-i]=-(etabinsA[mxBetaStrMod+i]=wMK->mSmdEGeom->EtaB()[i]);
86  for(int i=0;i<mxBMod2Pi+1;i++)
87  etaphibinsA[i]=wMK->mSmdEGeom->PhiB()[i];
88  hBsmdEtaAdc=new TH2F("eveBsmdEtaAdc"," Event: BSMD-Eta ADC vs. eta & phi; pseudorapidity; azimuth",mxBetaStrMod*2,etabinsA,mxBMod2Pi,etaphibinsA);
89 #endif
90 
91 }
92 
93 
94 
95 //-----------------------------
96 //-----------------------------
97 void
98 WeventDisplay::clear(){
99  hEmcET->Reset();
100  hTpcET->Reset();
101  for(int iep=0;iep<mxBSmd;iep++) hBsmdAdc[iep]->Reset();
102  for(int iuv=0;iuv<mxEsmdPlane;iuv++) hEsmdShower[iuv]->Reset();
103  hEsmdXpt->Reset();
104 
105 }
106 
107 //-----------------------------
108 //-----------------------------
109 void
110 WeventDisplay::draw( const char *tit,int eveID, int daqSeq, int runNo, WeveVertex myV, WeveEleTrack myTr){
111  if(maxEve<=0) return;
112  maxEve--;
113  TStyle *myStyle=new TStyle();
114  myStyle->cd();
115  myStyle->SetPalette(1,0);
116  myStyle->SetOptStat(1000010);
117 
118  char txt[1000];
119  sprintf(txt,"display-%s_run%d.eventId%06dvert%d",tit,runNo,eveID,myV.id);
120 
121  printf("WeventDisplay::Draw %s\n",txt);
122  TCanvas *c0; TPaveText *pvt;
123  string detector=tit;
124 
125  if(detector.compare("WB") == 0){ //barrel event display
126  sprintf(txt,"display-%s%.0f_run%d.eventId%05dvert%d",tit,myTr.cluster.ET,runNo,eveID,myV.id);
127  TFile hf(Form("%s.root",txt),"recreate"); //TFile hf(txt,"recreate");
128  c0=new TCanvas(txt,txt,850,600);
129  c0->cd();
130 
131  TString tt=txt;
132  TPad *cU = new TPad(tt+"U", tt+"U",0.,0.2,1.,1.); cU->Draw();
133  TPad *cD = new TPad(tt+"D", tt+"D",0.,0.,1.,0.2); cD->Draw();
134  cU->cd();
135  TPad *cU1 = new TPad(tt+"U1", tt+"U1",0.,0.,0.24,1.); cU1->Draw();
136  TPad *cU2 = new TPad(tt+"U2", tt+"U2",0.24,0.,0.55,1.); cU2->Draw();
137  TPad *cU3 = new TPad(tt+"U3", tt+"U3",0.55,0.,1.,1.); cU3->Draw();
138  cU3->Divide(2,1);
139  cU1->cd(); hTpcET->Draw("colz");
140 
141  TVector3 rW=myTr.pointTower.R;
142  rW.SetZ(rW.z()-myV.z);
143 
144  TEllipse *te1=new TEllipse(rW.Eta(),rW.Phi(), 0.1,0.1);
145  te1->SetFillStyle(0);te1->SetLineStyle(3); te1->SetLineColor(kMagenta);
146  TEllipse *te2=new TEllipse(rW.Eta(),rW.Phi(), 0.7, 0.7);
147  te2->SetFillStyle(0);te2->SetLineStyle(3); te2->SetLineColor(kBlack);
148 
149  TVector3 rA=-rW; // away direction
150  bxT->SetY1(rA.Phi() - wMK->par_awayDeltaPhi);
151  bxT->SetY2(rA.Phi() + wMK->par_awayDeltaPhi);
152 
153  bxE->SetY1(rA.Phi() - wMK->par_awayDeltaPhi);
154  bxE->SetY2(rA.Phi() + wMK->par_awayDeltaPhi);
155 
156 
157  te1->Draw(); te2->Draw(); bxT->Draw("l");
158  etaBL_ln->Draw(); etaBR_ln->Draw(); etaEL_ln->Draw();
159 
160  cU2->cd(); hEmcET->Draw("colz");
161  te1->Draw(); te2->Draw(); bxE->Draw("l");
162  etaBL_ln->Draw(); etaBR_ln->Draw();
163  etaEL_ln->Draw(); etaER_ln->Draw();
164 
165  for(int iep=0;iep<mxBSmd;iep++) {
166  cU3->cd(1+iep);
167  hBsmdAdc[iep]->Draw("colz");
168  te1->Draw(); te2->Draw(); bxT->Draw("l");
169  etaBL_ln->Draw(); etaBR_ln->Draw();
170  }
171 
172  //........... text information .............
173  pvt = new TPaveText(0,0.,1,1,"br");
174  cD->cd();
175  sprintf(txt,"run=%d eveID=%05d daq=%d vertex:ID=%d Z=%.0fcm",runNo,eveID,daqSeq,myV.id,myV.z);
176  printf("WeventDisplay::Event ID %s\n",txt);
177  pvt->AddText(txt);
178 
179  sprintf(txt,"TPC PT(GeV/c) near=%.1f away=%.1f ", myTr.nearTpcPT, myTr.awayTpcPT);
180  printf("WeventDisplay::Event TPC %s\n",txt);
181  pvt->AddText(txt);
182 
183  sprintf(txt,"BTOW ET/GeV: 2x2=%.1f near= %.1f away= %.1f",myTr.cluster.ET,myTr.nearBtowET,myTr.awayBtowET);
184  printf("WeventDisplay:: BTOW %s\n",txt);
185  pvt->AddText(txt);
186 
187  sprintf(txt,"Emc (Btow+Etow) ET/GeV: near= %.1f away= %.1f",myTr.nearEmcET,myTr.awayEmcET);
188  printf("WeventDisplay:: BTOW+ETOW %s\n",txt);
189  pvt->AddText(txt);
190 
191  sprintf(txt,"total ET/GeV: near= %.1f away= %.1f ptBalance= %.1f",myTr.nearTotET,myTr.awayTotET,myTr.ptBalance.Perp());
192  printf("WeventDisplay:: BTOW %s\n",txt);
193  pvt->AddText(txt);
194  // save dump of histos
195  if(hf.IsOpen()) {
196  hEmcET->Write();
197  hTpcET->Write();
198  for(int iep=0;iep<mxBSmd;iep++) hBsmdAdc[iep]->Write();
199  hf.Close();
200  }
201  }
202  else if(detector.compare("WE") == 0){ //endcap event display
203  sprintf(txt,"display-%s%.0f_run%d.eventId%05dvert%d",tit,myTr.cluster.ET,runNo,eveID,myV.id);
204  TFile hf(Form("%s.root",txt),"recreate");
205  c0=new TCanvas(txt,txt,1750,1300);
206  c0->cd();
207 
208  TString tt=txt;
209  TPad *cL = new TPad(tt+"L", tt+"L",0.,0.,0.6,1.); cL->Draw();
210  TPad *cR = new TPad(tt+"R", tt+"R",0.6,0.,1.,1.); cR->Draw();
211  cL->cd();
212  TPad *cLU = new TPad(tt+"LU", tt+"LU",0.,0.2,1.,1.); cLU->Draw();
213  TPad *cLD = new TPad(tt+"LD", tt+"LD",0.,0.,1.,0.2); cLD->Draw();
214  cLU->cd();
215  TPad *cLU1 = new TPad(tt+"LU1", tt+"LU1",0.,0.,0.44,1.); cLU1->Draw();
216  TPad *cLU2 = new TPad(tt+"LU2", tt+"LU2",0.44,0.,1.,1.); cLU2->Draw();
217  cR->cd();
218  TPad *cRU = new TPad(tt+"RU", tt+"RU",0.,0.5,1.,1.); cRU->Draw();
219  TPad *cRD = new TPad(tt+"RD", tt+"RD",0.,0.,1.,0.5); cRD->Draw();
220  cRD->Divide(1,2);
221  cLU1->cd(); hTpcET->Draw("colz");
222 
223  TVector3 rW=myTr.pointTower.R;
224  rW.SetZ(rW.z()-myV.z);
225  TEllipse *te1=new TEllipse(rW.Eta(),rW.Phi(), 0.1,0.1);
226  te1->SetFillStyle(0);te1->SetLineStyle(3); te1->SetLineColor(kMagenta);
227  TEllipse *te2=new TEllipse(rW.Eta(),rW.Phi(), 0.7, 0.7);
228  te2->SetFillStyle(0);te2->SetLineStyle(3); te2->SetLineColor(kBlack);
229 
230  TVector3 rA=-rW; // away direction
231  bxT->SetY1(rA.Phi() - wMK->par_awayDeltaPhi);
232  bxT->SetY2(rA.Phi() + wMK->par_awayDeltaPhi);
233 
234  bxE->SetY1(rA.Phi() - wMK->par_awayDeltaPhi);
235  bxE->SetY2(rA.Phi() + wMK->par_awayDeltaPhi);
236 
237  te1->Draw(); te2->Draw(); bxT->Draw("l");
238  etaBL_ln->Draw(); etaBR_ln->Draw(); etaEL_ln->Draw();
239 
240  cLU2->cd(); hEmcET->Draw("colz");
241  te1->Draw(); te2->Draw(); bxE->Draw("l");
242  etaBL_ln->Draw(); etaBR_ln->Draw();
243  etaEL_ln->Draw(); etaER_ln->Draw();
244 
245  TList *Lx; TLine *tline; TLine *tlineGlob;
246  //Draw ESMD showers
247  for(int iuv=0; iuv<mxEsmdPlane; iuv++){
248  cRD->cd(1+iuv);
249 
250  Lx=hEsmdShower[iuv]->GetListOfFunctions();
251  tline=new TLine(myTr.esmdDca[iuv],0.,myTr.esmdDca[iuv],100.);
252  tline->SetLineColor(2);
253  tlineGlob=new TLine(myTr.esmdGlobStrip[iuv]*0.5+myTr.esmdDcaGlob[iuv],0,myTr.esmdGlobStrip[iuv]*0.5+myTr.esmdDcaGlob[iuv],100.); tlineGlob->SetLineColor(2);
254  hEsmdShower[iuv]->Draw();
255  tline->Draw(); Lx->Add(tline);
256  tlineGlob->SetLineStyle(2); tlineGlob->Draw(); Lx->Add(tlineGlob);
257  hEsmdShower[iuv]->Draw();
258 
259  bxEs[iuv]->SetX1( 0.5*(-wMK->parE_esmdGL+myTr.esmdPeakOffset[iuv]) -0.25);
260  bxEs[iuv]->SetX2( 0.5*(+wMK->parE_esmdGL+myTr.esmdPeakOffset[iuv]+1)-0.25);
261  bxEs[iuv]->SetY2( myTr.esmdPeakSumE[iuv]/(2*wMK->parE_esmdGL+1));
262  bxEs[iuv]->Draw();
263 
264  //print Q/Pt warning
265  if( iuv==1 && myTr.prMuTrack->pt() >= 100.0 ) {
266  TLatex* tx = new TLatex(3,0.9*hEsmdShower[iuv]->GetMaximum(),"| Q/P_{T} | < 0.01"); tx->SetTextColor(kRed); tx->SetTextSize(0.1); tx->Draw();
267  }
268  }
269  hEsmdShower[0]->SetFillColor(kBlue-5);
270  hEsmdShower[1]->SetFillColor(kGreen-5);
271 
272  //make SMD cluster picture
273  cRU->cd();
274  hEsmdXpt->Draw("colz");
275  //fill with lines of hit strips
276  int secLoop[3]={myTr.hitSector-2,myTr.hitSector-1,myTr.hitSector};
277  if(secLoop[0] < 0) secLoop[0]=11;
278  if(secLoop[2] > 11) secLoop[2]=0;
279  for(int iuv=0; iuv<mxEsmdPlane; iuv++){
280  for(int isec=0; isec<3; isec++){
281  for(int k=0; k<288; k++){ //loop all strips
282  const StructEEmcStrip *stripPtr = wMK->geoSmd->getStripPtr(k,iuv,secLoop[isec]); //myTr.esmdStripId[iuv][k]
283  if(!stripPtr) continue;
284  if(wMK->wEve->esmd.ene[secLoop[isec]][iuv][k]*1e3 < 10.) continue;
285  int nColor = (int) (wMK->wEve->esmd.ene[secLoop[isec]][iuv][k]*1e3)/20;
286  int nSub=-10;
287  if(nColor==0) nSub=10;
288  if(nColor==1) nSub=8;
289  if(nColor==2) nSub=5;
290  if(nColor==3) nSub=1;
291  if(nColor==4) nSub=-4;
292  TVector3 end1=stripPtr->end1;
293  TVector3 end2=stripPtr->end2;
294  Lx=hEsmdXpt->GetListOfFunctions();
295  tline=new TLine(end1.X(),end1.Y(),end2.X(),end2.Y());
296  if(iuv==0) tline->SetLineColor(kBlue-nSub);
297  if(iuv==1) tline->SetLineColor(kGreen-nSub);
298  if(wMK->wEve->esmd.adc[secLoop[isec]][iuv][k] > 3896) tline->SetLineColor(2); //saturation with ped=200
299  tline->Draw(); Lx->Add(tline);
300  }
301  }
302  }
303 
304  TEllipse *te3=new TEllipse(rW.X(),rW.Y(), 2.0, 2.0);
305  te3->SetLineColor(kBlue); te3->SetFillStyle(0); te3->Draw(); //track projection
306  //TEllipse *te5=new TEllipse(myTr.pointTower.Rglob.X(),myTr.pointTower.Rglob.Y(), .25, .25);
307  //te5->SetLineColor(kBlue); te5->Draw(); //global track projection
308  TEllipse *te4=new TEllipse(0.,0., 215., 215.);
309  te4->SetFillStyle(0); te4->SetLineColor(kBlue); te4->Draw("l"); //outer esmd radius
310  TEllipse *te6=new TEllipse(0.,0., 75., 75.);
311  te6->SetFillStyle(0); te6->SetLineColor(kBlue); te6->Draw("l"); //inner esmd radius
312 
313  for(int iphi=0; iphi<12; iphi++){
314  float phi=(15.+ iphi*30.)/180.*TMath::Pi();
315  TVector3 r; r.SetPtThetaPhi(215,0,phi);
316  TVector3 rIn; rIn.SetPtThetaPhi(75,0,phi);
317  TLine *tl=new TLine(rIn.X(),rIn.Y(),r.X(),r.Y());
318  tl->SetLineColor(kBlue); tl->Draw(); //sector boundaries
319  }
320 
321  //........... text information .............
322  pvt = new TPaveText(0,0.,1,1,"br");
323  TH1F* hText = new TH1F("text"," ",1,0,1);
324  cLD->cd();
325  sprintf(txt,"run=%d eveID=%05d daq=%d vertex:ID=%d Z=%.0fcm ",runNo,eveID,daqSeq,myV.id,myV.z);
326  printf("WeventDisplay::Event ID %s\n",txt);
327  pvt->AddText(txt); hText->SetTitle(Form("%s%s",hText->GetTitle(),txt));
328  sprintf(txt,"TPC PT(GeV/c) prim=%.1f near=%.1f away=%.1f ", myTr.prMuTrack->pt(),myTr.nearTpcPT, myTr.awayTpcPT);
329  printf("WeventDisplay::Event %s\n",txt);
330  pvt->AddText(txt); hText->SetTitle(Form("%s%s",hText->GetTitle(),txt));
331 
332  sprintf(txt,"ETOW ET/GeV: 2x2=%.1f EMC: near= %.1f away= %.1f ",myTr.cluster.ET,myTr.nearEmcET,myTr.awayEmcET);
333  printf("WeventDisplay:: %s\n",txt);
334  pvt->AddText(txt); hText->SetTitle(Form("%s%s",hText->GetTitle(),txt));
335 
336  sprintf(txt,"total ET/GeV: near= %.1f away= %.1f ptBalance= %.1f ",myTr.nearTotET,myTr.awayTotET,myTr.ptBalance.Perp());
337  printf("WeventDisplay:: %s\n",txt);
338  pvt->AddText(txt); hText->SetTitle(Form("%s%s",hText->GetTitle(),txt));
339 
340  sprintf(txt,"Q/Pt = %.3f : ESMD E/MeV U peak= %.1f V peak= %.1f ",(1.0*myTr.prMuTrack->charge())/myTr.prMuTrack->pt(),myTr.esmdPeakSumE[0],myTr.esmdPeakSumE[1]);
341  printf("WeventDisplay:: %s\n",txt);
342  pvt->AddText(txt); hText->SetTitle(Form("%s%s",hText->GetTitle(),txt));
343 
344  float chi2=myTr.glMuTrack->chi2(); if(chi2>999.) chi2=-1.;
345  sprintf(txt,"Track: eta=%.1f Q=%d nFit=%d nPoss=%d r1=%.0f r2=%.0f chi2=%.1f",myTr.pointTower.R.Eta(),myTr.prMuTrack->charge(),myTr.prMuTrack->nHitsFit(),myTr.prMuTrack->nHitsPoss(),myTr.glMuTrack->firstPoint().perp(),myTr.glMuTrack->lastPoint().perp(),chi2);
346  printf("WeventDisplay:: %s\n",txt);
347  pvt->AddText(txt); hText->SetTitle(Form("%s%s",hText->GetTitle(),txt));
348 
349  // save dump of histos
350  if(hf.IsOpen()) {
351  hText->Write();
352  hEmcET->Write();
353  hTpcET->Write();
354  hEsmdShower[0]->Write();
355  hEsmdShower[1]->Write();
356  hEsmdXpt->Write();
357  hf.Close();
358  }
359  }
360 
361  pvt->Draw();
362  c0->Print();
363 
364 }
365 
366 //-----------------------------
367 //-----------------------------
368 void
369 WeventDisplay::exportEvent( const char *tit, WeveVertex myV, WeveEleTrack myTr, int vertexIndex){
370  if(maxEve<=0) return;
371  clear();
372  int eveId=wMK->wEve->id; //wMK->mMuDstMaker->muDst()->event()->eventId();
373  int runNo=wMK->wEve->runNo; //wMK->mMuDstMaker->muDst()->event()->runId();
374  const char *afile = ""; //wMK->mMuDstMaker->GetFile();
375  int len=strlen(afile);
376  int daqSeq=atoi(afile+(len-18));
377  // printf("DDD %s len=%d %d =%s=\n",afile,len,daqSeq,afile+(len-15));
378 
379  TVector3 rTw=myTr.cluster.position;
380  rTw.SetZ(rTw.z()-myV.z);
381  //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,
382  // runNo,daqSeq,eveId,myV.id,myV.z,
383  // myTr.prMuTrack->id(),myTr.prMuTrack->eta(),myTr.prMuTrack->phi()/3.1416*180.,myTr.glMuTrack->pt(),
384  // myTr.pointTower.id,wMK->wEve->bemc.adcTile[kBTow][myTr.pointTower.id-1],
385  // rTw.Eta(),rTw.Phi()/3.1416*180.,myTr.cluster.ET);
386 
387  float zVert=myV.z;
388  printf("WeventDisplay-%s::export run=%d eve=%d\n",tit,runNo,eveId);
389 
390  //.... process BTOW hits
391  for(int i=0;i< mxBtow;i++) {
392  float ene=wMK->wEve->bemc.eneTile[kBTow][i];
393  if(ene<=0) continue;
394  TVector3 primP=wMK->positionBtow[i]-TVector3(0,0,zVert);
395  primP.SetMag(ene); // it is 3D momentum in the event ref frame
396  float ET=primP.Perp();
397 
398  float eveEta=primP.Eta();
399  float evePhi=primP.Phi();
400  hEmcET->Fill(eveEta,evePhi,ET);
401  }
402 
403  //.... store ETOW hits
404  for(int i=0;i<mxEtowPhiBin;i++){
405  for(int j=0;j<mxEtowEta;j++){
406  float ene=wMK->wEve->etow.ene[i][j];
407  if(ene<=0) continue;
408  TVector3 primP=wMK->positionEtow[i][j]-TVector3(0,0,zVert);
409  primP.SetMag(ene); // it is 3D momentum in the event ref frame
410  float ET=primP.Perp();
411 
412  float eveEta=primP.Eta();
413  float evePhi=primP.Phi();
414  hEmcET->Fill(eveEta,evePhi,ET);
415  }
416  }
417 
418  hEmcET->SetMinimum(0.3); hEmcET->SetMaximum(30.);
419  // compute approximate event eta for barrel
420  float x,y,z;
421  float Rcylinder= wMK->mBtowGeom->Radius();
422  assert(wMK->mBtowGeom->getXYZ(20,x,y,z)==0); // this is approximate Z of last tower
423  TVector3 rL(Rcylinder,0,z+myV.z);
424  TVector3 rR(Rcylinder,0,z-myV.z);
425  float etaL=-rL.Eta(), etaR=rR.Eta();
426  etaBL_ln->SetX1(etaL); etaBL_ln->SetX2(etaL);
427  etaBR_ln->SetX1(etaR); etaBR_ln->SetX2(etaR);
428 
429  // compute approximate event eta for Endcap
430  rL=TVector3(0,214,270-myV.z); // detector eta~1.06
431  rR=TVector3(0,77,270-myV.z); // detector eta 2.0
432  etaL=rL.Eta(); etaR=rR.Eta();
433  etaEL_ln->SetX1(etaL); etaEL_ln->SetX2(etaL);
434  etaER_ln->SetX1(etaR); etaER_ln->SetX2(etaR);
435 
436 
437  //... TPC
438  hTpcET->SetMinimum(0.3);hTpcET->SetMaximum(10.);
439  if(wMK->mMuDstMaker) getPrimTracks( myV.id,myTr.pointTower.id);
440  else getPrimTracksFromTree(vertexIndex,myTr.pointTower.id);
441 
442  //.... BSMD-Eta, -Phi
443 
444  for(int iep=0;iep<mxBSmd;iep++) {
445  hBsmdAdc[iep]->SetMinimum(30);hBsmdAdc[iep]->SetMaximum(999);
446  for(int i=0;i< mxBStrips;i++) {
447  float adc=wMK->wEve->bemc.adcBsmd[iep][i];
448  if(adc<=0) continue;
449  TVector3 r=wMK->positionBsmd[iep][i];
450  float z1=r.z()-zVert;
451  r.SetZ(z1);
452  hBsmdAdc[iep]->Fill(r.Eta(),r.Phi(),adc);
453  }
454  }// end of eta,phi-planes
455 
456  //.... ESMD shower shape
457  for(int iuv=0; iuv<mxEsmdPlane; iuv++)
458  for(int j=0;j<41;j++)
459  hEsmdShower[iuv]->SetBinContent(j+1,myTr.esmdShower[iuv][j]);
460 
461  //... ESMD cluster picture
462  //initialize histo centred at track extrapolation
463  TVector3 rW=myTr.pointTower.R; //z is at SMD depth
464  float width=65.;
465  hEsmdXpt->SetBins(130,rW.X()-width,rW.X()+width,130,rW.Y()-width,rW.Y()+width);
466 
467  //.... produce plot & save
468  draw(tit,eveId, daqSeq,runNo,myV, myTr);
469  //export2sketchup(tit,myV, myTr);
470 }
471 
472 //-----------------------------
473 //-----------------------------
474 void
475 WeventDisplay::getPrimTracks( int vertID,int pointTowId) {
476  assert(vertID>=0);
477  assert(vertID<(int)wMK->mMuDstMaker->muDst()->numberOfPrimaryVertices());
478  StMuPrimaryVertex* V=wMK-> mMuDstMaker->muDst()->primaryVertex(vertID);
479  assert(V);
480  wMK-> mMuDstMaker->muDst()->setVertexIndex(vertID);
481  float rank=V->ranking();
482  assert(rank>0);
483 
484  Int_t nPrimTrAll=wMK->mMuDstMaker->muDst()->GetNPrimaryTrack();
485  for(int itr=0;itr<nPrimTrAll;itr++) {
486  StMuTrack *prTr=wMK->mMuDstMaker->muDst()->primaryTracks(itr);
487  if(prTr->flag()<=0) continue;
488  if(prTr->flag()!=301 && pointTowId>0) continue;// TPC-only regular tracks for barrel candidate
489  if(prTr->flag()!=301 && prTr->flag()!=311 && pointTowId<0) continue;// TPC regular and short EEMC tracks for endcap candidate
490  float hitFrac=1.*prTr->nHitsFit()/prTr->nHitsPoss();
491  if(hitFrac<wMK->par_nHitFrac) continue;
492  StThreeVectorF prPvect=prTr->p();
493  TVector3 primP=TVector3(prPvect.x(),prPvect.y(),prPvect.z());
494  float pT=prTr->pt();
495  hTpcET->Fill(prTr->eta(),prTr->phi(),pT);
496 
497  }
498 }
499 
500 
501 //-----------------------------
502 //-----------------------------
503 void
504 WeventDisplay::getPrimTracksFromTree(int vertID,int pointTowId) {
505 
506  // flag=2 use 2D cut, 1= only delta phi
507 
508  assert(vertID>=0);
509  assert(vertID<(int)wMK->wEve->vertex.size());
510 
511  WeveVertex &V=wMK->wEve->vertex[vertID];
512  for(unsigned int it=0;it<V.prTrList.size();it++){
513  StMuTrack *prTr=V.prTrList[it];
514  if(prTr->flag()<=0) continue;
515  if(prTr->flag()!=301 && pointTowId>0) continue;// TPC-only regular tracks for barrel candidate
516  if(prTr->flag()!=301 && prTr->flag()!=311 && pointTowId<0) continue;// TPC regular and short EEMC tracks for endcap candidate
517  float hitFrac=1.*prTr->nHitsFit()/prTr->nHitsPoss();
518  if(hitFrac<wMK->par_nHitFrac) continue;
519  StThreeVectorF prPvect=prTr->p();
520  TVector3 primP=TVector3(prPvect.x(),prPvect.y(),prPvect.z());
521  float pT=prTr->pt();
522  hTpcET->Fill(prTr->eta(),prTr->phi(),pT);
523 
524  }
525 }
526 
527 
528 //-----------------------------
529 //-----------------------------
530 void
531 WeventDisplay::export2sketchup( const char *tit, WeveVertex myV, WeveEleTrack myTr){
532  int eveId=wMK->mMuDstMaker->muDst()->event()->eventId();
533  int runNo=wMK->mMuDstMaker->muDst()->event()->runId();
534  char txt[1000];
535  sprintf(txt,"display3D-%s_run%d.eventId%05d_vert%d.txt",tit,runNo,eveId,myV.id);
536  FILE *fd=fopen(txt,"w"); assert(fd);
537 
538  //........ DUMP PRIM TRACKS..........
539  int vertID=myV.id;
540  assert(vertID>=0);
541  assert(vertID<(int)wMK->mMuDstMaker->muDst()->numberOfPrimaryVertices());
542  StMuPrimaryVertex* V=wMK-> mMuDstMaker->muDst()->primaryVertex(vertID);
543  assert(V);
544  wMK-> mMuDstMaker->muDst()->setVertexIndex(vertID);
545  float rank=V->ranking();
546  assert(rank>0);
547  const StThreeVectorF &rV=V->position();
548  Int_t nPrimTrAll=wMK->mMuDstMaker->muDst()->GetNPrimaryTrack();
549  for(int itr=0;itr<nPrimTrAll;itr++) {
550  StMuTrack *prTr=wMK->mMuDstMaker->muDst()->primaryTracks(itr);
551  if(prTr->flag()<=0) continue;
552  if(prTr->flag()!=301) continue;// TPC-only regular tracks
553  float hitFrac=1.*prTr->nHitsFit()/prTr->nHitsPoss();
554  if(hitFrac<wMK->par_nHitFrac) continue;
555  prTr->p();
556  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());
557  }
558 
559  //........DUMP BTOW towers
560  float Rcylinder= wMK->mBtowGeom->Radius(), Rcylinder2=Rcylinder*Rcylinder;
561  for(int i=0;i< mxBtow;i++) {
562  float ene=wMK->wEve->bemc.eneTile[kBTow][i];
563  if(ene<=0) continue;
564  float delZ=wMK->positionBtow[i].z()-myV.z;
565  float e2et=Rcylinder/sqrt(Rcylinder2+delZ*delZ);
566  float ET=ene*e2et;
567  float detEta=wMK->positionBtow[i].Eta();
568  float detPhi=wMK->positionBtow[i].Phi();
569  fprintf(fd,"btow V %.1f %.3f %.3f eveET:detEta:detPhi %.3f %.3f %.3f\n",rV.x(),rV.y(),rV.z(),ET, detEta,detPhi);
570  }
571 
572  const char cPlane[ mxBSmd]={'E','P'};
573  //........DUMP BSMD hits
574  for(int iep=0;iep<mxBSmd;iep++) {
575  for(int i=0;i< mxBStrips;i++) {
576  float adc=wMK->wEve->bemc.adcBsmd[iep][i];
577  if(adc<=0) continue;
578  TVector3 r=wMK->positionBsmd[iep][i];
579  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() );
580  }
581  }
582 
583  //........DUMP ETOW towers
584  for(int iphi=0; iphi<mxEtowPhiBin; iphi++){
585  for(int ieta=0; ieta<mxEtowEta; ieta++){//sum all eta rings
586  float ene=wMK->wEve->etow.ene[iphi][ieta];
587  if(ene<=0) continue; //skip towers with no energy
588  TVector3 detP=wMK->positionEtow[iphi][ieta];
589  TVector3 primP=detP-TVector3(0,0,myV.z);
590  primP.SetMag(ene); // it is 3D momentum in the event ref frame
591  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());
592  }
593  }
594 
595  //.... DUMP reco electron ...
596  float eleET=myTr.cluster.ET;
597  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());
598 
599 
600  //... close event file
601  fclose(fd);
602 
603 
604 }
605 
606 
607 // $Log: WeventDisplay.cxx,v $
608 // Revision 1.6 2012/09/21 16:59:10 balewski
609 // added ESMD peak adjustement - partialy finished
610 //
611 // Revision 1.5 2012/09/18 21:10:09 stevens4
612 // Include all rank>0 vertex again (new jet format coming next), and remove rank<0 endcap vertices.
613 //
614 // Revision 1.4 2012/08/21 18:29:16 stevens4
615 // Updates to endcap W selection using ESMD strip ratio
616 //
617 // Revision 1.3 2012/06/25 20:53:29 stevens4
618 // algo and histo cleanup
619 //
620 // Revision 1.2 2012/06/18 18:28:01 stevens4
621 // Updates for Run 9+11+12 AL analysis
622 //
623 // Revision 1.1 2011/02/10 20:33:26 balewski
624 // start
625 //
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::chi2
Double_t chi2() const
Returns chi2 of fit.
Definition: StMuTrack.h:251
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
StMuTrack::firstPoint
const StThreeVectorF & firstPoint() const
Returns positions of first measured point.
Definition: StMuTrack.h:261
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::lastPoint
const StThreeVectorF & lastPoint() const
Returns positions of last measured point.
Definition: StMuTrack.h:262
StMuTrack::phi
Double_t phi() const
Returns phi at point of dca to primary vertex.
Definition: StMuTrack.h:258
St2011WMaker
muDst based extraction of W-signal from pp500 data from 2011
Definition: St2011WMaker.h:49
EEmcSmdGeom::getStripPtr
StructEEmcStrip * getStripPtr(const Int_t iStrip, const Int_t iUV, const Int_t iSec)
return a strip pointer from indices
Definition: EEmcSmdGeom.cxx:375
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