StRoot: StSpinPool/StWalgoB2009/WeventDisplay.cxx Source File

00001 // $Id: WeventDisplay.cxx,v 1.6 2010/01/10 03:01:37 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 #include "St2009WMaker.h"
00030 #include "WanaConst.h"
00031 #include "WeventDisplay.h"
00032 //-----------------------------
00033 //-----------------------------
00034 WeventDisplay::WeventDisplay( St2009WMaker* mk, int mxEv) {
00035   maxEve=mxEv;
00036   wMK=mk;
00037   const float PI=TMath::Pi();
00038   const char cPlane[ mxBSmd]={'E','P'};
00039   char txt1[100], txt2[1000];
00040 
00041   // barrel
00042   etaBL_ln=new TLine(-1,-3.2,1,3.2); etaBL_ln->SetLineColor(kBlue);
00043   etaBR_ln=new TLine(-1,-3.2,1,3.2); etaBR_ln->SetLineColor(kBlue);
00044 
00045   etaEL_ln=new TLine(2,-3.2,3,3.2); etaEL_ln->SetLineColor(kGreen);
00046   etaER_ln=new TLine(2,-3.2,3,3.2); etaER_ln->SetLineColor(kGreen);
00047 
00048   bxT=new TBox(-1.3,0, 1.3,1.);  bxT->SetFillStyle(0); bxT->SetLineStyle(2);
00049   bxE=new TBox(-1.3,0, 1,2.);  bxE->SetFillStyle(0); bxE->SetLineStyle(2);
00050   if(wMK->par_useEtow>=2) bxE->SetX2(2.2);
00051 
00052 
00053   hEmcET=new TH2F("eveBtowET","EMC  ET sum, Z=[0.3,30]GeV; event eta ; phi",38,-1.4,2.4,63,-PI,PI);
00054 
00055   hTpcET=new TH2F("eveTpcET","TPC PT sum, Z[0.3,10]GeV/c; event eta ; phi",28,-1.4,1.4,63,-PI,PI);
00056 
00057   for(int iep=0;iep<mxBSmd;iep++) {
00058     sprintf(txt1,"eveBsmdAdc_%c",cPlane[iep]);
00059     sprintf(txt2,"BSMD_%c  ADC sum Z=[30,1000]; event eta ; phi",cPlane[iep]);
00060     hBsmdAdc[iep]=new TH2F(txt1,txt2,26,-1.3,1.3,63,-PI,PI);
00061   }
00062 
00063   
00064 #if 0
00065   //---- smart code from Willie
00066   float etabinsA[1+mxBetaStrMod*2], etaphibinsA[mxBMod2Pi+1];
00067   for(int i=0;i<mxBetaStrMod+1;i++) 
00068     etabinsA[mxBetaStrMod-i]=-(etabinsA[mxBetaStrMod+i]=wMK->mSmdEGeom->EtaB()[i]);
00069   for(int i=0;i<mxBMod2Pi+1;i++) 
00070     etaphibinsA[i]=wMK->mSmdEGeom->PhiB()[i];
00071   hBsmdEtaAdc=new TH2F("eveBsmdEtaAdc"," Event: BSMD-Eta ADC vs. eta & phi; pseudorapidity; azimuth",mxBetaStrMod*2,etabinsA,mxBMod2Pi,etaphibinsA);
00072 #endif  
00073   
00074 }
00075 
00076 
00077 
00078 //-----------------------------
00079 //-----------------------------
00080 void
00081 WeventDisplay::clear(){
00082   hEmcET->Reset();
00083   hTpcET->Reset();
00084   for(int iep=0;iep<mxBSmd;iep++)  hBsmdAdc[iep]->Reset();
00085 
00086 }
00087 
00088 //-----------------------------
00089 //-----------------------------
00090 void
00091 WeventDisplay::draw(  const char *tit,int eveID, int daqSeq,  int runNo,  WeveVertex myV, WeveEleTrack myTr){
00092   if(maxEve<=0) return;
00093   maxEve--;
00094   TStyle *myStyle=new TStyle();
00095   myStyle->cd();
00096   myStyle->SetPalette(1,0);
00097   myStyle->SetOptStat(1000010);
00098   
00099   char txt[1000];
00100   sprintf(txt,"display-%s_run%d.eventId%06dvert%d",tit,runNo,eveID,myV.id);
00101  
00102   printf("WeventDisplay::Draw %s\n",txt);
00103   TCanvas *c0=new TCanvas(txt,txt,850,600);
00104   c0->cd();
00105 
00106   TString tt=txt;
00107   TPad *cU = new TPad(tt+"U", tt+"U",0.,0.2,1.,1.); cU->Draw();
00108   TPad *cD = new TPad(tt+"D", tt+"D",0.,0.,1.,0.2); cD->Draw();
00109 
00110   cU->cd();
00111   TPad *cU1 = new TPad(tt+"U1", tt+"U1",0.,0.,0.24,1.); cU1->Draw();
00112   TPad *cU2 = new TPad(tt+"U2", tt+"U2",0.24,0.,0.55,1.); cU2->Draw();
00113   TPad *cU3 = new TPad(tt+"U3", tt+"U3",0.55,0.,1.,1.);  cU3->Draw();
00114   cU3->Divide(2,1);
00115  
00116 
00117   cU1->cd(); hTpcET->Draw("colz");
00118      
00119   TVector3 rW=myTr.pointTower.R;
00120   rW.SetZ(rW.z()-myV.z);
00121 
00122   TEllipse *te1=new TEllipse(rW.Eta(),rW.Phi(), 0.1,0.1);
00123   te1->SetFillStyle(0);te1->SetLineStyle(3); te1->SetLineColor(kMagenta);
00124   TEllipse *te2=new TEllipse(rW.Eta(),rW.Phi(), 0.7, 0.7);
00125   te2->SetFillStyle(0);te2->SetLineStyle(3); te2->SetLineColor(kBlack);
00126 
00127   TVector3 rA=-rW; // away direction
00128   bxT->SetY1(rA.Phi() - wMK->par_awayDeltaPhi);
00129   bxT->SetY2(rA.Phi() + wMK->par_awayDeltaPhi);
00130 
00131   bxE->SetY1(rA.Phi() - wMK->par_awayDeltaPhi);
00132   bxE->SetY2(rA.Phi() + wMK->par_awayDeltaPhi);
00133 
00134 
00135   te1->Draw();   te2->Draw(); bxT->Draw("l");
00136   etaBL_ln->Draw();  etaBR_ln->Draw(); etaEL_ln->Draw();
00137   
00138   cU2->cd();    hEmcET->Draw("colz");
00139   te1->Draw();  te2->Draw(); bxE->Draw("l");
00140   etaBL_ln->Draw();  etaBR_ln->Draw();
00141   etaEL_ln->Draw();  etaER_ln->Draw();
00142   
00143   for(int iep=0;iep<mxBSmd;iep++) {
00144     cU3->cd(1+iep);
00145     hBsmdAdc[iep]->Draw("colz");
00146     te1->Draw();  te2->Draw();  bxT->Draw("l");
00147     etaBL_ln->Draw();  etaBR_ln->Draw();
00148   }
00149 
00150   //........... text information .............
00151   TPaveText *pvt = new TPaveText(0,0.,1,1,"br");
00152   cD->cd();
00153   sprintf(txt,"run=%d  eveID=%05d daq=%d vertex:ID=%d Z=%.0fcm",runNo,eveID,daqSeq,myV.id,myV.z);
00154   printf("WeventDisplay::Event ID  %s\n",txt);
00155   pvt->AddText(txt);
00156 
00157   sprintf(txt,"TPC PT(GeV/c) near=%.1f  away=%.1f ", myTr.nearTpcPT, myTr.awayTpcPT);
00158   printf("WeventDisplay::Event TPC  %s\n",txt);
00159   pvt->AddText(txt);
00160 
00161   sprintf(txt,"BTOW ET/GeV: 2x2=%.1f   near= %.1f   away= %.1f",myTr.cluster.ET,myTr.nearBtowET,myTr.awayBtowET);
00162   printf("WeventDisplay:: BTOW  %s\n",txt);
00163   pvt->AddText(txt);
00164 
00165   sprintf(txt,"Emc (Btow+Etow) ET/GeV:   near= %.1f   away= %.1f",myTr.nearEmcET,myTr.awayEmcET);
00166   printf("WeventDisplay:: BTOW+ETOW  %s\n",txt);
00167   pvt->AddText(txt);
00168 
00169   sprintf(txt,"total ET/GeV:   near= %.1f   away= %.1f  ptBalance= %.1f",myTr.nearTotET,myTr.awayTotET,myTr.ptBalance.Perp());
00170   printf("WeventDisplay:: BTOW  %s\n",txt);
00171   pvt->AddText(txt);
00172 
00173 
00174   pvt->Draw();
00175 
00176   c0->Print();
00177   // save dump of histos
00178   sprintf(txt,"display-%s_run%d.eventId%05dvert%d.root",tit,runNo,eveID,myV.id);
00179   TFile hf(txt,"recreate");
00180   if(hf.IsOpen()) {
00181     hEmcET->Write();
00182     hTpcET->Write();
00183     for(int iep=0;iep<mxBSmd;iep++) hBsmdAdc[iep]->Write();
00184     hf.Close();
00185   }
00186 }
00187 
00188 //-----------------------------
00189 //-----------------------------
00190 void
00191 WeventDisplay::exportEvent( const char *tit, WeveVertex myV, WeveEleTrack myTr){
00192   if(maxEve<=0) return;
00193   clear();
00194   int eveId=wMK->mMuDstMaker->muDst()->event()->eventId();
00195   int runNo=wMK->mMuDstMaker->muDst()->event()->runId();
00196   const char *afile = wMK->mMuDstMaker->GetFile();
00197   int len=strlen(afile);
00198   int daqSeq=atoi(afile+(len-18));
00199   //  printf("DDD %s len=%d %d =%s=\n",afile,len,daqSeq,afile+(len-15));
00200 
00201   TVector3 rTw=myTr.cluster.position;
00202   rTw.SetZ(rTw.z()-myV.z);
00203   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,
00204          runNo,daqSeq,eveId,myV.id,myV.z,
00205          myTr.prMuTrack->id(),myTr.prMuTrack->eta(),myTr.prMuTrack->phi()/3.1416*180.,myTr.glMuTrack->pt(),
00206          myTr.pointTower.id,wMK->wEve.bemc.adcTile[kBTow][myTr.pointTower.id-1],
00207          rTw.Eta(),rTw.Phi()/3.1416*180.,myTr.cluster.ET);
00208 
00209   float zVert=myV.z;
00210   printf("WeventDisplay-%s::export run=%d eve=%d\n",tit,runNo,eveId);
00211 
00212  //.... process BTOW hits
00213   for(int i=0;i< mxBtow;i++) {
00214     float ene=wMK->wEve.bemc.eneTile[kBTow][i];
00215     if(ene<=0) continue;
00216     TVector3 primP=wMK->positionBtow[i]-TVector3(0,0,zVert);
00217     primP.SetMag(ene); // it is 3D momentum in the event ref frame
00218     float ET=primP.Perp();
00219 
00220     float eveEta=primP.Eta();
00221     float evePhi=primP.Phi();
00222     hEmcET->Fill(eveEta,evePhi,ET); 
00223   }
00224   
00225   int par_useEtow = wMK->par_useEtow; //get ETOW flag   
00226   //.... store ETOW hits    
00227   if(par_useEtow >= 1) {
00228     for(int i=0;i<mxEtowPhiBin;i++){
00229       for(int j=0;j<mxEtowEta;j++){
00230         float ene=wMK->wEve.etow.ene[i][j];
00231         if(ene<=0) continue;
00232         TVector3 primP=wMK->positionEtow[i][j]-TVector3(0,0,zVert);
00233         primP.SetMag(ene); // it is 3D momentum in the event ref frame
00234         float ET=primP.Perp();
00235         
00236         float eveEta=primP.Eta();
00237         float evePhi=primP.Phi();
00238         hEmcET->Fill(eveEta,evePhi,ET); 
00239       }
00240     }
00241   }
00242 
00243   hEmcET->SetMinimum(0.3);  hEmcET->SetMaximum(30.);
00244  // compute approximate event eta for barrel
00245   float x,y,z;
00246   float Rcylinder= wMK->mBtowGeom->Radius();
00247   assert(wMK->mBtowGeom->getXYZ(20,x,y,z)==0); // this is approximate Z of last tower
00248   TVector3 rL(Rcylinder,0,z+myV.z);
00249   TVector3 rR(Rcylinder,0,z-myV.z);
00250   float etaL=-rL.Eta(), etaR=rR.Eta();
00251   etaBL_ln->SetX1(etaL);  etaBL_ln->SetX2(etaL);
00252   etaBR_ln->SetX1(etaR);  etaBR_ln->SetX2(etaR);
00253   if(wMK->par_useEtow<2) bxE->SetX2(etaR+0.05); // skip endcap
00254   
00255   // compute approximate event eta for Endcap
00256   rL=TVector3(0,214,270-myV.z); // detector eta~1.06
00257   rR=TVector3(0,77,270-myV.z); // detector eta 2.0
00258   etaL=rL.Eta(); etaR=rR.Eta();
00259   etaEL_ln->SetX1(etaL);  etaEL_ln->SetX2(etaL);
00260   etaER_ln->SetX1(etaR);  etaER_ln->SetX2(etaR);
00261 
00262 
00263   //... TPC
00264   hTpcET->SetMinimum(0.3);hTpcET->SetMaximum(10.);
00265   getPrimTracks( myV.id);
00266   
00267   //.... BSMD-Eta, -Phi
00268   
00269   for(int iep=0;iep<mxBSmd;iep++) {
00270     hBsmdAdc[iep]->SetMinimum(30);hBsmdAdc[iep]->SetMaximum(999);
00271     for(int i=0;i< mxBStrips;i++) {
00272       float adc=wMK->wEve.bemc.adcBsmd[iep][i];
00273       if(adc<=0) continue;
00274       TVector3 r=wMK->positionBsmd[iep][i];
00275       float z1=r.z()-zVert;
00276       r.SetZ(z1);
00277       hBsmdAdc[iep]->Fill(r.Eta(),r.Phi(),adc);
00278     }
00279   }// end of eta,phi-planes
00280 
00281   //.... produce plot & save
00282   draw(tit,eveId, daqSeq,runNo,myV, myTr);
00283   export2sketchup(tit,myV, myTr);
00284 }
00285 
00286 //-----------------------------
00287 //-----------------------------
00288 void
00289 WeventDisplay::getPrimTracks( int vertID) {
00290   assert(vertID>=0);
00291   assert(vertID<(int)wMK->mMuDstMaker->muDst()->numberOfPrimaryVertices());
00292   StMuPrimaryVertex* V=wMK-> mMuDstMaker->muDst()->primaryVertex(vertID);
00293   assert(V);
00294   wMK-> mMuDstMaker->muDst()->setVertexIndex(vertID);
00295   float rank=V->ranking();
00296   assert(rank>0);
00297   
00298   Int_t nPrimTrAll=wMK->mMuDstMaker->muDst()->GetNPrimaryTrack();
00299   for(int itr=0;itr<nPrimTrAll;itr++) {
00300     StMuTrack *prTr=wMK->mMuDstMaker->muDst()->primaryTracks(itr);
00301     if(prTr->flag()<=0) continue;
00302     if(prTr->flag()!=301) continue;// TPC-only regular tracks
00303     float hitFrac=1.*prTr->nHitsFit()/prTr->nHitsPoss();
00304     if(hitFrac<wMK->par_nHitFrac) continue;
00305     StThreeVectorF prPvect=prTr->p();
00306     TVector3 primP=TVector3(prPvect.x(),prPvect.y(),prPvect.z());
00307     float pT=prTr->pt();
00308     hTpcET->Fill(prTr->eta(),prTr->phi(),pT);
00309 
00310   }
00311 }
00312 
00313 
00314 //-----------------------------
00315 //-----------------------------
00316 void
00317 WeventDisplay::export2sketchup(  const char *tit, WeveVertex myV, WeveEleTrack myTr){
00318   int eveId=wMK->mMuDstMaker->muDst()->event()->eventId();
00319   int runNo=wMK->mMuDstMaker->muDst()->event()->runId();
00320   char txt[1000];
00321   sprintf(txt,"display3D-%s_run%d.eventId%05d_vert%d.txt",tit,runNo,eveId,myV.id);
00322   FILE *fd=fopen(txt,"w"); assert(fd);
00323   
00324   //........ DUMP PRIM TRACKS..........
00325   int vertID=myV.id;
00326   assert(vertID>=0);
00327   assert(vertID<(int)wMK->mMuDstMaker->muDst()->numberOfPrimaryVertices());
00328   StMuPrimaryVertex* V=wMK-> mMuDstMaker->muDst()->primaryVertex(vertID);
00329   assert(V);
00330   wMK-> mMuDstMaker->muDst()->setVertexIndex(vertID);
00331   float rank=V->ranking();
00332   assert(rank>0);
00333   const StThreeVectorF &rV=V->position();
00334   Int_t nPrimTrAll=wMK->mMuDstMaker->muDst()->GetNPrimaryTrack();
00335   for(int itr=0;itr<nPrimTrAll;itr++) {
00336     StMuTrack *prTr=wMK->mMuDstMaker->muDst()->primaryTracks(itr);
00337     if(prTr->flag()<=0) continue;
00338     if(prTr->flag()!=301) continue;// TPC-only regular tracks
00339     float hitFrac=1.*prTr->nHitsFit()/prTr->nHitsPoss();
00340     if(hitFrac<wMK->par_nHitFrac) continue;
00341     prTr->p();
00342     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());
00343   }
00344 
00345   //........DUMP BTOW towers
00346   float Rcylinder= wMK->mBtowGeom->Radius(), Rcylinder2=Rcylinder*Rcylinder;
00347   for(int i=0;i< mxBtow;i++) {
00348     float ene=wMK->wEve.bemc.eneTile[kBTow][i];
00349     if(ene<=0) continue;
00350     float delZ=wMK->positionBtow[i].z()-myV.z;
00351     float e2et=Rcylinder/sqrt(Rcylinder2+delZ*delZ);
00352     float ET=ene*e2et;    
00353     float detEta=wMK->positionBtow[i].Eta();
00354     float detPhi=wMK->positionBtow[i].Phi();
00355     fprintf(fd,"btow V %.1f %.3f %.3f  eveET:detEta:detPhi %.3f %.3f  %.3f\n",rV.x(),rV.y(),rV.z(),ET, detEta,detPhi);
00356   }
00357   
00358   const char cPlane[ mxBSmd]={'E','P'};
00359   //........DUMP BSMD  hits 
00360   for(int iep=0;iep<mxBSmd;iep++) {
00361     for(int i=0;i< mxBStrips;i++) {
00362       float adc=wMK->wEve.bemc.adcBsmd[iep][i];
00363       if(adc<=0) continue;
00364       TVector3 r=wMK->positionBsmd[iep][i];
00365       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() );
00366     }
00367   }
00368 
00369   //........DUMP ETOW towers
00370   for(int iphi=0; iphi<mxEtowPhiBin; iphi++){
00371     for(int ieta=0; ieta<mxEtowEta; ieta++){//sum all eta rings
00372       float ene=wMK->wEve.etow.ene[iphi][ieta];
00373       if(ene<=0) continue; //skip towers with no energy
00374       TVector3 detP=wMK->positionEtow[iphi][ieta];
00375       TVector3 primP=detP-TVector3(0,0,myV.z);
00376       primP.SetMag(ene); // it is 3D momentum in the event ref frame
00377       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());
00378     }
00379   }
00380   
00381   //.... DUMP reco electron ...
00382   float eleET=myTr.cluster.ET;
00383   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()); 
00384 
00385   
00386   //... close event file
00387   fclose(fd);
00388 
00389 
00390 }
00391 
00392 
00393 // $Log: WeventDisplay.cxx,v $
00394 // Revision 1.6  2010/01/10 03:01:37  balewski
00395 // cleanup & nicer histos
00396 //
00397 // Revision 1.5  2010/01/09 00:07:16  stevens4
00398 // add jet finder
00399 //
00400 // Revision 1.4  2010/01/06 19:16:48  stevens4
00401 // track cuts now on primary component, cleanup
00402 //
00403 // Revision 1.3  2010/01/06 04:22:15  balewski
00404 // added Q/PT plot for Zs, more cleanup
00405 //
00406 // Revision 1.2  2009/12/10 16:01:31  stevens4
00407 // fixed zero-length vector
00408 //
00409 // Revision 1.1  2009/11/23 23:00:18  balewski
00410 // code moved spin-pool
00411 //