StRoot: StSpinPool/StWalgo2011/WeventDisplay.cxx Source File

00001 // $Id: WeventDisplay.cxx,v 1.1 2011/02/10 20:33:26 balewski Exp $
00002 //
00003 //*-- Author : Jan Balewski, MIT
00004 
00005 
00006 #include <TH1.h>
00007 #include <TH2.h>
00008 #include <TMath.h>
00009 #include <TCanvas.h>
00010 #include <TStyle.h>
00011 #include <TFile.h>
00012 #include <TText.h>
00013 #include <TList.h>
00014 #include <TBox.h>
00015 #include <TPaveText.h>
00016 #include <TPad.h>
00017 #include <TEllipse.h>
00018 #include <stdio.h>
00019 
00020 //MuDst
00021 #include <StMuDSTMaker/COMMON/StMuDstMaker.h>
00022 #include <StMuDSTMaker/COMMON/StMuDst.h>
00023 #include <StMuDSTMaker/COMMON/StMuTriggerIdCollection.h>
00024 #include <StMuDSTMaker/COMMON/StMuEvent.h>
00025 #include <StMuDSTMaker/COMMON/StMuTrack.h>
00026 #include <StMuDSTMaker/COMMON/StMuPrimaryVertex.h>
00027 #include "StEmcUtil/geometry/StEmcGeom.h"
00028 
00029 //Esmd geometry for plots
00030 #include "StEEmcUtil/EEmcGeom/EEmcGeomSimple.h"
00031 #include "StEEmcUtil/StEEmcSmd/EEmcSmdGeom.h"
00032 
00033 #include "St2011WMaker.h"
00034 #include "WanaConst.h"
00035 #include "WeventDisplay.h"
00036 //-----------------------------
00037 //-----------------------------
00038 WeventDisplay::WeventDisplay( St2011WMaker* mk, int mxEv) {
00039   maxEve=mxEv;
00040   wMK=mk;
00041   const float PI=TMath::Pi();
00042   const char cPlane[ mxBSmd]={'E','P'};
00043   const char cEsmdPlane[ mxEsmdPlane]={'U','V'};
00044   char txt1[100], txt2[1000];
00045 
00046   // barrel
00047   etaBL_ln=new TLine(-1,-3.2,1,3.2); etaBL_ln->SetLineColor(kBlue);
00048   etaBR_ln=new TLine(-1,-3.2,1,3.2); etaBR_ln->SetLineColor(kBlue);
00049 
00050   etaEL_ln=new TLine(2,-3.2,3,3.2); etaEL_ln->SetLineColor(kGreen);
00051   etaER_ln=new TLine(2,-3.2,3,3.2); etaER_ln->SetLineColor(kGreen);
00052 
00053   bxT=new TBox(-1.3,0, 1.3,1.);  bxT->SetFillStyle(0); bxT->SetLineStyle(2);
00054   bxE=new TBox(-1.3,0, 1,2.);  bxE->SetFillStyle(0); bxE->SetLineStyle(2);
00055   bxE->SetX2(2.2);
00056 
00057 
00058   hEmcET=new TH2F("eveBtowET","EMC  ET sum, Z=[0.3,30]GeV; event eta ; phi",38,-1.4,2.4,63,-PI,PI);
00059 
00060   hTpcET=new TH2F("eveTpcET","TPC PT sum, Z[0.3,10]GeV/c; event eta ; phi",28,-1.4,1.4,63,-PI,PI);
00061 
00062   for(int iep=0;iep<mxBSmd;iep++) {
00063     sprintf(txt1,"eveBsmdAdc_%c",cPlane[iep]);
00064     sprintf(txt2,"BSMD_%c  ADC sum Z=[30,1000]; event eta ; phi",cPlane[iep]);
00065     hBsmdAdc[iep]=new TH2F(txt1,txt2,26,-1.3,1.3,63,-PI,PI);
00066   }
00067 
00068   //ESMD shower shape
00069   for(int iuv=0; iuv<mxEsmdPlane; iuv++){
00070     sprintf(txt1,"eveEsmdShower_%c",cEsmdPlane[iuv]);
00071     sprintf(txt2,"ESMD_%c Shower Shape; i_strip position - track position (cm) ; MeV",cEsmdPlane[iuv]);
00072     hEsmdShower[iuv]=new TH1F(txt1,txt2,41,-10.25,10.25);
00073   }
00074   hEsmdXpt=new TH2F("eveEsmdXpt","ESMD Cross Point; X (cm); Y (cm)",100,0.,100,100,0.,100.);
00075 
00076   
00077 #if 0
00078   //---- smart code from Willie
00079   float etabinsA[1+mxBetaStrMod*2], etaphibinsA[mxBMod2Pi+1];
00080   for(int i=0;i<mxBetaStrMod+1;i++) 
00081     etabinsA[mxBetaStrMod-i]=-(etabinsA[mxBetaStrMod+i]=wMK->mSmdEGeom->EtaB()[i]);
00082   for(int i=0;i<mxBMod2Pi+1;i++) 
00083     etaphibinsA[i]=wMK->mSmdEGeom->PhiB()[i];
00084   hBsmdEtaAdc=new TH2F("eveBsmdEtaAdc"," Event: BSMD-Eta ADC vs. eta & phi; pseudorapidity; azimuth",mxBetaStrMod*2,etabinsA,mxBMod2Pi,etaphibinsA);
00085 #endif  
00086   
00087 }
00088 
00089 
00090 
00091 //-----------------------------
00092 //-----------------------------
00093 void
00094 WeventDisplay::clear(){
00095   hEmcET->Reset();
00096   hTpcET->Reset();
00097   for(int iep=0;iep<mxBSmd;iep++)  hBsmdAdc[iep]->Reset();
00098   for(int iuv=0;iuv<mxEsmdPlane;iuv++)  hEsmdShower[iuv]->Reset();
00099   hEsmdXpt->Reset();
00100 
00101 }
00102 
00103 //-----------------------------
00104 //-----------------------------
00105 void
00106 WeventDisplay::draw(  const char *tit,int eveID, int daqSeq,  int runNo,  WeveVertex myV, WeveEleTrack myTr){
00107   if(maxEve<=0) return;
00108   maxEve--;
00109   TStyle *myStyle=new TStyle();
00110   myStyle->cd();
00111   myStyle->SetPalette(1,0);
00112   myStyle->SetOptStat(1000010);
00113   
00114   char txt[1000];
00115   sprintf(txt,"display-%s_run%d.eventId%06dvert%d",tit,runNo,eveID,myV.id);
00116  
00117   printf("WeventDisplay::Draw %s\n",txt);
00118   TCanvas *c0; TPaveText *pvt;
00119   string detector=tit;
00120 
00121   if(detector.compare("WB") == 0){ //barrel event display
00122     c0=new TCanvas(txt,txt,850,600);
00123     c0->cd();
00124     
00125     TString tt=txt;
00126     TPad *cU = new TPad(tt+"U", tt+"U",0.,0.2,1.,1.); cU->Draw();
00127     TPad *cD = new TPad(tt+"D", tt+"D",0.,0.,1.,0.2); cD->Draw();
00128     cU->cd();
00129     TPad *cU1 = new TPad(tt+"U1", tt+"U1",0.,0.,0.24,1.); cU1->Draw();
00130     TPad *cU2 = new TPad(tt+"U2", tt+"U2",0.24,0.,0.55,1.); cU2->Draw();
00131     TPad *cU3 = new TPad(tt+"U3", tt+"U3",0.55,0.,1.,1.);  cU3->Draw();
00132     cU3->Divide(2,1);
00133     cU1->cd(); hTpcET->Draw("colz");
00134     
00135     TVector3 rW=myTr.pointTower.R;
00136     rW.SetZ(rW.z()-myV.z);
00137     
00138     TEllipse *te1=new TEllipse(rW.Eta(),rW.Phi(), 0.1,0.1);
00139     te1->SetFillStyle(0);te1->SetLineStyle(3); te1->SetLineColor(kMagenta);
00140     TEllipse *te2=new TEllipse(rW.Eta(),rW.Phi(), 0.7, 0.7);
00141     te2->SetFillStyle(0);te2->SetLineStyle(3); te2->SetLineColor(kBlack);
00142     
00143     TVector3 rA=-rW; // away direction
00144     bxT->SetY1(rA.Phi() - wMK->par_awayDeltaPhi);
00145     bxT->SetY2(rA.Phi() + wMK->par_awayDeltaPhi);
00146     
00147     bxE->SetY1(rA.Phi() - wMK->par_awayDeltaPhi);
00148     bxE->SetY2(rA.Phi() + wMK->par_awayDeltaPhi);
00149     
00150     
00151     te1->Draw();   te2->Draw(); bxT->Draw("l");
00152     etaBL_ln->Draw();  etaBR_ln->Draw(); etaEL_ln->Draw();
00153     
00154     cU2->cd();    hEmcET->Draw("colz");
00155     te1->Draw();  te2->Draw(); bxE->Draw("l");
00156     etaBL_ln->Draw();  etaBR_ln->Draw();
00157     etaEL_ln->Draw();  etaER_ln->Draw();
00158     
00159     for(int iep=0;iep<mxBSmd;iep++) {
00160       cU3->cd(1+iep);
00161       hBsmdAdc[iep]->Draw("colz");
00162       te1->Draw();  te2->Draw();  bxT->Draw("l");
00163       etaBL_ln->Draw();  etaBR_ln->Draw();
00164     }
00165     
00166     //........... text information .............
00167     pvt = new TPaveText(0,0.,1,1,"br");
00168     cD->cd();
00169     sprintf(txt,"run=%d  eveID=%05d daq=%d vertex:ID=%d Z=%.0fcm",runNo,eveID,daqSeq,myV.id,myV.z);
00170     printf("WeventDisplay::Event ID  %s\n",txt);
00171     pvt->AddText(txt);
00172     
00173     sprintf(txt,"TPC PT(GeV/c) near=%.1f  away=%.1f ", myTr.nearTpcPT, myTr.awayTpcPT);
00174     printf("WeventDisplay::Event TPC  %s\n",txt);
00175     pvt->AddText(txt);
00176     
00177     sprintf(txt,"BTOW ET/GeV: 2x2=%.1f   near= %.1f   away= %.1f",myTr.cluster.ET,myTr.nearBtowET,myTr.awayBtowET);
00178     printf("WeventDisplay:: BTOW  %s\n",txt);
00179     pvt->AddText(txt);
00180     
00181     sprintf(txt,"Emc (Btow+Etow) ET/GeV:   near= %.1f   away= %.1f",myTr.nearEmcET,myTr.awayEmcET);
00182     printf("WeventDisplay:: BTOW+ETOW  %s\n",txt);
00183     pvt->AddText(txt);
00184     
00185     sprintf(txt,"total ET/GeV:   near= %.1f   away= %.1f  ptBalance= %.1f",myTr.nearTotET,myTr.awayTotET,myTr.ptBalance.Perp());
00186     printf("WeventDisplay:: BTOW  %s\n",txt);
00187     pvt->AddText(txt);
00188     // save dump of histos
00189     sprintf(txt,"display-%s_run%d.eventId%05dvert%d.root",tit,runNo,eveID,myV.id);
00190     TFile hf(txt,"recreate");
00191     if(hf.IsOpen()) {
00192       hEmcET->Write();
00193       hTpcET->Write();
00194       for(int iep=0;iep<mxBSmd;iep++) hBsmdAdc[iep]->Write();
00195       hf.Close();
00196     }
00197   }
00198   else if(detector.compare("WE") == 0){ //endcap event display
00199     c0=new TCanvas(txt,txt,1750,1300);
00200     c0->cd();
00201 
00202     TString tt=txt;
00203     TPad *cL = new TPad(tt+"L", tt+"L",0.,0.,0.6,1.); cL->Draw();
00204     TPad *cR = new TPad(tt+"R", tt+"R",0.6,0.,1.,1.); cR->Draw();
00205     cL->cd();
00206     TPad *cLU = new TPad(tt+"LU", tt+"LU",0.,0.2,1.,1.); cLU->Draw();
00207     TPad *cLD = new TPad(tt+"LD", tt+"LD",0.,0.,1.,0.2); cLD->Draw();
00208     cLU->cd();
00209     TPad *cLU1 = new TPad(tt+"LU1", tt+"LU1",0.,0.,0.44,1.); cLU1->Draw();
00210     TPad *cLU2 = new TPad(tt+"LU2", tt+"LU2",0.44,0.,1.,1.); cLU2->Draw();
00211     cR->cd();
00212     TPad *cRU = new TPad(tt+"RU", tt+"RU",0.,0.5,1.,1.); cRU->Draw();
00213     TPad *cRD = new TPad(tt+"RD", tt+"RD",0.,0.,1.,0.5); cRD->Draw();
00214     cRD->Divide(1,2);
00215     cLU1->cd(); hTpcET->Draw("colz");
00216     
00217     TVector3 rW=myTr.pointTower.R;
00218     rW.SetZ(rW.z()-myV.z);
00219     TEllipse *te1=new TEllipse(rW.Eta(),rW.Phi(), 0.1,0.1);
00220     te1->SetFillStyle(0);te1->SetLineStyle(3); te1->SetLineColor(kMagenta);
00221     TEllipse *te2=new TEllipse(rW.Eta(),rW.Phi(), 0.7, 0.7);
00222     te2->SetFillStyle(0);te2->SetLineStyle(3); te2->SetLineColor(kBlack);
00223     
00224     TVector3 rA=-rW; // away direction
00225     bxT->SetY1(rA.Phi() - wMK->par_awayDeltaPhi);
00226     bxT->SetY2(rA.Phi() + wMK->par_awayDeltaPhi);
00227     
00228     bxE->SetY1(rA.Phi() - wMK->par_awayDeltaPhi);
00229     bxE->SetY2(rA.Phi() + wMK->par_awayDeltaPhi);
00230     
00231     te1->Draw();   te2->Draw(); bxT->Draw("l");
00232     etaBL_ln->Draw();  etaBR_ln->Draw(); etaEL_ln->Draw();
00233     
00234     cLU2->cd();    hEmcET->Draw("colz");
00235     te1->Draw();  te2->Draw(); bxE->Draw("l");
00236     etaBL_ln->Draw();  etaBR_ln->Draw();
00237     etaEL_ln->Draw();  etaER_ln->Draw();
00238     
00239     TList *Lx;  TLine *tline; TLine *tlineGlob;
00240     //Draw ESMD showers
00241     for(int iuv=0; iuv<mxEsmdPlane; iuv++){
00242       cRD->cd(1+iuv);
00243 
00244       Lx=hEsmdShower[iuv]->GetListOfFunctions();
00245       tline=new TLine(myTr.esmdDca[iuv],0.,myTr.esmdDca[iuv],100.);
00246       tline->SetLineColor(2);
00247       tlineGlob=new TLine(myTr.esmdGlobStrip[iuv]*0.5+myTr.esmdDcaGlob[iuv],0,myTr.esmdGlobStrip[iuv]*0.5+myTr.esmdDcaGlob[iuv],100.); tlineGlob->SetLineColor(2);
00248       hEsmdShower[iuv]->Draw();
00249       tline->Draw(); Lx->Add(tline);
00250       tlineGlob->SetLineStyle(2); tlineGlob->Draw(); Lx->Add(tlineGlob);
00251       hEsmdShower[iuv]->Draw();
00252     }
00253     hEsmdShower[0]->SetFillColor(kBlue-5);
00254     hEsmdShower[1]->SetFillColor(kGreen-5);
00255 
00256     //make SMD cluster picture
00257     cRU->cd();
00258     hEsmdXpt->Draw("colz");
00259     //fill with lines of hit strips
00260     int secLoop[3]={myTr.hitSector-2,myTr.hitSector-1,myTr.hitSector};
00261     if(secLoop[0] < 0) secLoop[0]=11;
00262     if(secLoop[2] > 11) secLoop[2]=0;
00263     for(int iuv=0; iuv<mxEsmdPlane; iuv++){
00264       for(int isec=0; isec<3; isec++){
00265         for(int k=0; k<288; k++){ //loop all strips
00266           const StructEEmcStrip *stripPtr = wMK->geoSmd->getStripPtr(k,iuv,secLoop[isec]); //myTr.esmdStripId[iuv][k]
00267           if(!stripPtr) continue;
00268           if(wMK->wEve->esmd.ene[secLoop[isec]][iuv][k]*1e3 < 10.) continue;
00269           int nColor = (int) (wMK->wEve->esmd.ene[secLoop[isec]][iuv][k]*1e3)/20;
00270           int nSub=-10;
00271           if(nColor==0) nSub=10;
00272           if(nColor==1) nSub=8;
00273           if(nColor==2) nSub=5;
00274           if(nColor==3) nSub=1;
00275           if(nColor==4) nSub=-4;
00276           TVector3 end1=stripPtr->end1;
00277           TVector3 end2=stripPtr->end2;
00278           Lx=hEsmdXpt->GetListOfFunctions();
00279           tline=new TLine(end1.X(),end1.Y(),end2.X(),end2.Y());
00280           if(iuv==0) tline->SetLineColor(kBlue-nSub);
00281           if(iuv==1) tline->SetLineColor(kGreen-nSub);
00282           if(wMK->wEve->esmd.ene[secLoop[isec]][iuv][k]*1e3 > 100) tline->SetLineColor(2);
00283           tline->Draw(); Lx->Add(tline);
00284         }
00285       }
00286     }
00287 
00288     TEllipse *te3=new TEllipse(rW.X(),rW.Y(), 2.0, 2.0);
00289     te3->SetLineColor(kBlue); te3->SetFillStyle(0); te3->Draw(); //track projection
00290     //TEllipse *te5=new TEllipse(myTr.pointTower.Rglob.X(),myTr.pointTower.Rglob.Y(), .25, .25);
00291     //te5->SetLineColor(kBlue); te5->Draw(); //global track projection
00292     TEllipse *te4=new TEllipse(0.,0., 215., 215.);
00293     te4->SetFillStyle(0); te4->SetLineColor(kBlue); te4->Draw("l"); //outer esmd radius
00294     TEllipse *te6=new TEllipse(0.,0., 75., 75.);
00295     te6->SetFillStyle(0); te6->SetLineColor(kBlue); te6->Draw("l"); //inner esmd radius
00296     
00297     for(int iphi=0; iphi<12; iphi++){
00298       float phi=(15.+ iphi*30.)/180.*TMath::Pi();
00299       TVector3 r; r.SetPtThetaPhi(215,0,phi);
00300       TVector3 rIn; rIn.SetPtThetaPhi(75,0,phi);
00301       TLine *tl=new TLine(rIn.X(),rIn.Y(),r.X(),r.Y());
00302       tl->SetLineColor(kBlue); tl->Draw(); //sector boundaries
00303     }
00304     
00305     //........... text information .............
00306     pvt = new TPaveText(0,0.,1,1,"br");
00307     cLD->cd();
00308     sprintf(txt,"run=%d  eveID=%05d daq=%d vertex:ID=%d Z=%.0fcm",runNo,eveID,daqSeq,myV.id,myV.z);
00309     printf("WeventDisplay::Event ID  %s\n",txt);
00310     pvt->AddText(txt);
00311     sprintf(txt,"TPC PT(GeV/c) near=%.1f  away=%.1f ", myTr.nearTpcPT, myTr.awayTpcPT);
00312     printf("WeventDisplay::Event %s\n",txt);
00313     pvt->AddText(txt);
00314     
00315     sprintf(txt,"ETOW ET/GeV: 2x1=%.1f  EMC: near= %.1f   away= %.1f",myTr.cluster.ET,myTr.nearEmcET,myTr.awayEmcET);
00316     printf("WeventDisplay:: %s\n",txt);
00317     pvt->AddText(txt);
00318     
00319     sprintf(txt,"total ET/GeV:   near= %.1f   away= %.1f  ptBalance= %.1f",myTr.nearTotET,myTr.awayTotET,myTr.ptBalance.Perp());
00320     printf("WeventDisplay:: %s\n",txt);
00321     pvt->AddText(txt);
00322     
00323     sprintf(txt,"ESMD E/MeV:   U plane= %.1f  V plane= %.1f ",myTr.esmdE[0],myTr.esmdE[1]);
00324     printf("WeventDisplay:: %s\n",txt);
00325     pvt->AddText(txt);
00326 
00327     float chi2=myTr.glMuTrack->chi2(); if(chi2>999.) chi2=-1.;
00328     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);
00329     printf("WeventDisplay:: %s\n",txt);
00330     pvt->AddText(txt);
00331   }
00332   
00333   pvt->Draw();
00334   c0->Print();
00335     
00336 }
00337 
00338 //-----------------------------
00339 //-----------------------------
00340 void
00341 WeventDisplay::exportEvent( const char *tit, WeveVertex myV, WeveEleTrack myTr){
00342   if(maxEve<=0) return;
00343   clear();
00344   int eveId=wMK->mMuDstMaker->muDst()->event()->eventId();
00345   int runNo=wMK->mMuDstMaker->muDst()->event()->runId();
00346   const char *afile = wMK->mMuDstMaker->GetFile();
00347   int len=strlen(afile);
00348   int daqSeq=atoi(afile+(len-18));
00349   //  printf("DDD %s len=%d %d =%s=\n",afile,len,daqSeq,afile+(len-15));
00350 
00351   TVector3 rTw=myTr.cluster.position;
00352   rTw.SetZ(rTw.z()-myV.z);
00353   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,
00354          runNo,daqSeq,eveId,myV.id,myV.z,
00355          myTr.prMuTrack->id(),myTr.prMuTrack->eta(),myTr.prMuTrack->phi()/3.1416*180.,myTr.glMuTrack->pt(),
00356          myTr.pointTower.id,wMK->wEve->bemc.adcTile[kBTow][myTr.pointTower.id-1],
00357          rTw.Eta(),rTw.Phi()/3.1416*180.,myTr.cluster.ET);
00358 
00359   float zVert=myV.z;
00360   printf("WeventDisplay-%s::export run=%d eve=%d\n",tit,runNo,eveId);
00361 
00362  //.... process BTOW hits
00363   for(int i=0;i< mxBtow;i++) {
00364     float ene=wMK->wEve->bemc.eneTile[kBTow][i];
00365     if(ene<=0) continue;
00366     TVector3 primP=wMK->positionBtow[i]-TVector3(0,0,zVert);
00367     primP.SetMag(ene); // it is 3D momentum in the event ref frame
00368     float ET=primP.Perp();
00369 
00370     float eveEta=primP.Eta();
00371     float evePhi=primP.Phi();
00372     hEmcET->Fill(eveEta,evePhi,ET); 
00373   }
00374   
00375   //.... store ETOW hits    
00376   for(int i=0;i<mxEtowPhiBin;i++){
00377     for(int j=0;j<mxEtowEta;j++){
00378       float ene=wMK->wEve->etow.ene[i][j];
00379       if(ene<=0) continue;
00380       TVector3 primP=wMK->positionEtow[i][j]-TVector3(0,0,zVert);
00381       primP.SetMag(ene); // it is 3D momentum in the event ref frame
00382       float ET=primP.Perp();
00383       
00384       float eveEta=primP.Eta();
00385       float evePhi=primP.Phi();
00386       hEmcET->Fill(eveEta,evePhi,ET);   
00387     }
00388   }
00389   
00390   hEmcET->SetMinimum(0.3);  hEmcET->SetMaximum(30.);
00391  // compute approximate event eta for barrel
00392   float x,y,z;
00393   float Rcylinder= wMK->mBtowGeom->Radius();
00394   assert(wMK->mBtowGeom->getXYZ(20,x,y,z)==0); // this is approximate Z of last tower
00395   TVector3 rL(Rcylinder,0,z+myV.z);
00396   TVector3 rR(Rcylinder,0,z-myV.z);
00397   float etaL=-rL.Eta(), etaR=rR.Eta();
00398   etaBL_ln->SetX1(etaL);  etaBL_ln->SetX2(etaL);
00399   etaBR_ln->SetX1(etaR);  etaBR_ln->SetX2(etaR);
00400     
00401   // compute approximate event eta for Endcap
00402   rL=TVector3(0,214,270-myV.z); // detector eta~1.06
00403   rR=TVector3(0,77,270-myV.z); // detector eta 2.0
00404   etaL=rL.Eta(); etaR=rR.Eta();
00405   etaEL_ln->SetX1(etaL);  etaEL_ln->SetX2(etaL);
00406   etaER_ln->SetX1(etaR);  etaER_ln->SetX2(etaR);
00407 
00408 
00409   //... TPC
00410   hTpcET->SetMinimum(0.3);hTpcET->SetMaximum(10.);
00411   getPrimTracks( myV.id);
00412   
00413   //.... BSMD-Eta, -Phi
00414   
00415   for(int iep=0;iep<mxBSmd;iep++) {
00416     hBsmdAdc[iep]->SetMinimum(30);hBsmdAdc[iep]->SetMaximum(999);
00417     for(int i=0;i< mxBStrips;i++) {
00418       float adc=wMK->wEve->bemc.adcBsmd[iep][i];
00419       if(adc<=0) continue;
00420       TVector3 r=wMK->positionBsmd[iep][i];
00421       float z1=r.z()-zVert;
00422       r.SetZ(z1);
00423       hBsmdAdc[iep]->Fill(r.Eta(),r.Phi(),adc);
00424     }
00425   }// end of eta,phi-planes
00426 
00427   //.... ESMD shower shape
00428   for(int iuv=0; iuv<mxEsmdPlane; iuv++)
00429     for(int j=0;j<41;j++)
00430       hEsmdShower[iuv]->SetBinContent(j+1,myTr.esmdShower[iuv][j]);
00431 
00432   //...  ESMD cluster picture
00433   //initialize histo centred at track extrapolation
00434   TVector3 rW=myTr.pointTower.R; //z is at SMD depth
00435   float width=65.;
00436   hEsmdXpt->SetBins(130,rW.X()-width,rW.X()+width,130,rW.Y()-width,rW.Y()+width);
00437 
00438   //.... produce plot & save
00439   draw(tit,eveId, daqSeq,runNo,myV, myTr);
00440   export2sketchup(tit,myV, myTr);
00441 }
00442 
00443 //-----------------------------
00444 //-----------------------------
00445 void
00446 WeventDisplay::getPrimTracks( int vertID) {
00447   assert(vertID>=0);
00448   assert(vertID<(int)wMK->mMuDstMaker->muDst()->numberOfPrimaryVertices());
00449   StMuPrimaryVertex* V=wMK-> mMuDstMaker->muDst()->primaryVertex(vertID);
00450   assert(V);
00451   wMK-> mMuDstMaker->muDst()->setVertexIndex(vertID);
00452   float rank=V->ranking();
00453   assert(rank>0 || (rank<0 && V->nEEMCMatch()));
00454   
00455   Int_t nPrimTrAll=wMK->mMuDstMaker->muDst()->GetNPrimaryTrack();
00456   for(int itr=0;itr<nPrimTrAll;itr++) {
00457     StMuTrack *prTr=wMK->mMuDstMaker->muDst()->primaryTracks(itr);
00458     if(prTr->flag()<=0) continue;
00459     if(prTr->flag()!=301) continue;// TPC-only regular tracks
00460     float hitFrac=1.*prTr->nHitsFit()/prTr->nHitsPoss();
00461     if(hitFrac<wMK->par_nHitFrac) continue;
00462     StThreeVectorF prPvect=prTr->p();
00463     TVector3 primP=TVector3(prPvect.x(),prPvect.y(),prPvect.z());
00464     float pT=prTr->pt();
00465     hTpcET->Fill(prTr->eta(),prTr->phi(),pT);
00466 
00467   }
00468 }
00469 
00470 
00471 //-----------------------------
00472 //-----------------------------
00473 void
00474 WeventDisplay::export2sketchup(  const char *tit, WeveVertex myV, WeveEleTrack myTr){
00475   int eveId=wMK->mMuDstMaker->muDst()->event()->eventId();
00476   int runNo=wMK->mMuDstMaker->muDst()->event()->runId();
00477   char txt[1000];
00478   sprintf(txt,"display3D-%s_run%d.eventId%05d_vert%d.txt",tit,runNo,eveId,myV.id);
00479   FILE *fd=fopen(txt,"w"); assert(fd);
00480   
00481   //........ DUMP PRIM TRACKS..........
00482   int vertID=myV.id;
00483   assert(vertID>=0);
00484   assert(vertID<(int)wMK->mMuDstMaker->muDst()->numberOfPrimaryVertices());
00485   StMuPrimaryVertex* V=wMK-> mMuDstMaker->muDst()->primaryVertex(vertID);
00486   assert(V);
00487   wMK-> mMuDstMaker->muDst()->setVertexIndex(vertID);
00488   float rank=V->ranking();
00489   assert(rank>0 || (rank<0 && V->nEEMCMatch()));
00490   const StThreeVectorF &rV=V->position();
00491   Int_t nPrimTrAll=wMK->mMuDstMaker->muDst()->GetNPrimaryTrack();
00492   for(int itr=0;itr<nPrimTrAll;itr++) {
00493     StMuTrack *prTr=wMK->mMuDstMaker->muDst()->primaryTracks(itr);
00494     if(prTr->flag()<=0) continue;
00495     if(prTr->flag()!=301) continue;// TPC-only regular tracks
00496     float hitFrac=1.*prTr->nHitsFit()/prTr->nHitsPoss();
00497     if(hitFrac<wMK->par_nHitFrac) continue;
00498     prTr->p();
00499     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());
00500   }
00501 
00502   //........DUMP BTOW towers
00503   float Rcylinder= wMK->mBtowGeom->Radius(), Rcylinder2=Rcylinder*Rcylinder;
00504   for(int i=0;i< mxBtow;i++) {
00505     float ene=wMK->wEve->bemc.eneTile[kBTow][i];
00506     if(ene<=0) continue;
00507     float delZ=wMK->positionBtow[i].z()-myV.z;
00508     float e2et=Rcylinder/sqrt(Rcylinder2+delZ*delZ);
00509     float ET=ene*e2et;    
00510     float detEta=wMK->positionBtow[i].Eta();
00511     float detPhi=wMK->positionBtow[i].Phi();
00512     fprintf(fd,"btow V %.1f %.3f %.3f  eveET:detEta:detPhi %.3f %.3f  %.3f\n",rV.x(),rV.y(),rV.z(),ET, detEta,detPhi);
00513   }
00514   
00515   const char cPlane[ mxBSmd]={'E','P'};
00516   //........DUMP BSMD  hits 
00517   for(int iep=0;iep<mxBSmd;iep++) {
00518     for(int i=0;i< mxBStrips;i++) {
00519       float adc=wMK->wEve->bemc.adcBsmd[iep][i];
00520       if(adc<=0) continue;
00521       TVector3 r=wMK->positionBsmd[iep][i];
00522       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() );
00523     }
00524   }
00525 
00526   //........DUMP ETOW towers
00527   for(int iphi=0; iphi<mxEtowPhiBin; iphi++){
00528     for(int ieta=0; ieta<mxEtowEta; ieta++){//sum all eta rings
00529       float ene=wMK->wEve->etow.ene[iphi][ieta];
00530       if(ene<=0) continue; //skip towers with no energy
00531       TVector3 detP=wMK->positionEtow[iphi][ieta];
00532       TVector3 primP=detP-TVector3(0,0,myV.z);
00533       primP.SetMag(ene); // it is 3D momentum in the event ref frame
00534       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());
00535     }
00536   }
00537   
00538   //.... DUMP reco electron ...
00539   float eleET=myTr.cluster.ET;
00540   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()); 
00541 
00542   
00543   //... close event file
00544   fclose(fd);
00545 
00546 
00547 }
00548 
00549 
00550 // $Log: WeventDisplay.cxx,v $
00551 // Revision 1.1  2011/02/10 20:33:26  balewski
00552 // start
00553 //