00001 #include "StClusterDisplay.h"
00002 #include "StEventTypes.h"
00003 #include "StEvent.h"
00004 #include "StTrack.h"
00005 #include <math.h>
00006 #include "TFile.h"
00007 #include "TROOT.h"
00008 #include "StEmcUtil/geometry/StEmcGeom.h"
00009 #include "StEmcUtil/others/emcDetectorName.h"
00010
00011 ClassImp(StClusterDisplay)
00012
00013
00014 StClusterDisplay::StClusterDisplay(const char *name):StMaker(name)
00015 {}
00016
00017 StClusterDisplay::~StClusterDisplay()
00018 {}
00019
00020 Int_t StClusterDisplay::Init()
00021 {
00022 TString name,title;
00023 for(Int_t i = 0; i<MAXDETBARREL; i++)
00024 {
00025 mGeo[i]=StEmcGeom::getEmcGeom(detname[i].Data());
00026 mCanvas[i] = 0;
00027 mTh[i] = 0.2;
00028
00029 name = detname[i];
00030 name += "_NClusters";
00031 title = "Number of clusters for detector ";
00032 title+= detname[i];
00033 mHist1D[0][i] = new TH1F(name.Data(),title.Data(),500,0,500);
00034
00035 name = detname[i];
00036 name += "_NHits";
00037 title = "Number of hits/cluster for detector ";
00038 title+= detname[i];
00039 mHist1D[1][i] = new TH1F(name.Data(),title.Data(),20,0,20);
00040
00041 name = detname[i];
00042 name += "_Energy";
00043 title = "Cluster energy for detector ";
00044 title+= detname[i];
00045 mHist1D[2][i] = new TH1F(name.Data(),title.Data(),500,0,50);
00046
00047 name = detname[i];
00048 name += "_RMSEta";
00049 title = "Cluster RMS/eta for detector ";
00050 title+= detname[i];
00051 mHist1D[3][i] = new TH1F(name.Data(),title.Data(),20,0,0.1);
00052
00053 name = detname[i];
00054 name += "_RMSPhi";
00055 title = "Cluster RMS/phi for detector ";
00056 title+= detname[i];
00057 mHist1D[4][i] = new TH1F(name.Data(),title.Data(),20,0,0.1);
00058
00059 name = detname[i];
00060 name += "_EtaPhi";
00061 title = "Cluster (eta,phi) for detector ";
00062 title+= detname[i];
00063 mHist2D[0][i] = new TH2F(name.Data(),title.Data(),40,-1,1,120,-3.1415,3.1415);
00064
00065
00066 Int_t nEta=mGeo[i]->NEta();
00067 Int_t nSub=mGeo[i]->NSub();
00068
00069 TArrayF EtaB(nEta+1,mGeo[i]->EtaB());
00070 TArrayF PhiB(nSub+1,mGeo[i]->PhiB());
00071
00072 TArrayF EtaBins(2*nEta+2);
00073 for(Int_t j=0;j<2*nEta+2;j++)
00074 {
00075 if (j<nEta+1)
00076 EtaBins[j]=-EtaB[nEta-j];
00077 else
00078 EtaBins[j]=EtaB[j-nEta-1];
00079 }
00080
00081 TArrayF PhiBins1(60*(nSub+1));
00082 TArrayF PhiBins(60*(nSub+1));
00083 Int_t j=0;
00084 for(Int_t m=1;m<=60;m++)
00085 for(Int_t s=1;s<=nSub+1;s++)
00086 {
00087 Float_t center;
00088 mGeo[i]->getPhiModule(m,center);
00089 PhiBins1[j]=center+PhiB[s-1];
00090 j++;
00091 }
00092
00093 Bool_t again=kTRUE;
00094 j=0;
00095 do
00096 {
00097 again=kFALSE;
00098 Float_t phitmp=6.4;
00099 Int_t ktmp=-1;
00100 for(Int_t k=0;k<60*(nSub+1);k++)
00101 {
00102 if(PhiBins1[k]<phitmp)
00103 {
00104 phitmp=PhiBins1[k];
00105 ktmp=k;
00106 }
00107 }
00108 if(ktmp!=-1)
00109 {
00110 PhiBins[j]=phitmp;
00111 again=kTRUE;
00112 PhiBins1[ktmp]=999;
00113 j++;
00114 }
00115 }
00116 while(again);
00117
00118 name = detname[i];
00119 name += "_HitDisplay";
00120 title = "Hits distribution for detector ";
00121 title+= detname[i];
00122 mHist2D[1][i] = new TH2F(name.Data(),title.Data(),2*nEta+2-1,EtaBins.GetArray(),60*(nSub+1)-1,PhiBins.GetArray());
00123
00124 name = detname[i];
00125 name += "_ClusterDisplay";
00126 title = "Cluster distribution for detector ";
00127 title+= detname[i];
00128 Int_t BF=10;
00129 mHist2D[2][i] = new TH2F(name.Data(),title.Data(),BF*2*nEta,-1,1,BF*60*nSub,-3.1415,3.1415);
00130 }
00131 return StMaker::Init();
00132 }
00133
00134 Int_t StClusterDisplay::Make()
00135 {
00136 StEvent *ev = (StEvent*)GetInputDS("StEvent");
00137 if(!ev)
00138 return kStOk;
00139
00140 StEmcCollection *emc = ev->emcCollection();
00141 if(!emc)
00142 return kStOk;
00143 for(Int_t i = 0;i<MAXDETBARREL;i++)
00144 {
00145 mHist2D[1][i]->Reset();
00146 mHist2D[2][i]->Reset();
00147 StDetectorId id = static_cast<StDetectorId>(i+kBarrelEmcTowerId);
00148 StEmcDetector* detector = emc->detector(id);
00149 if(detector)
00150 {
00151 for(UInt_t j=1;j<121;j++)
00152 {
00153 StEmcModule* module = detector->module(j);
00154 if(module)
00155 {
00156 StSPtrVecEmcRawHit& rawHit=module->hits();
00157
00158 for(Int_t k=0;k<(Int_t)rawHit.size();k++)
00159 {
00160 Float_t eta,phi;
00161 Int_t m=rawHit[k]->module();
00162 Int_t e=rawHit[k]->eta();
00163 Int_t s=abs(rawHit[k]->sub());
00164 Float_t energy = rawHit[k]->energy();
00165
00166 mGeo[i]->getEta(m,e,eta);
00167 mGeo[i]->getPhi(m,s,phi);
00168 if(energy>mTh[i])
00169 mHist2D[1][i]->Fill(eta,phi,energy);
00170 }
00171 }
00172 }
00173 StEmcClusterCollection* coll = detector->cluster();
00174 if(coll)
00175 {
00176 StSPtrVecEmcCluster& clusters = coll->clusters();
00177 Int_t n = clusters.size();
00178 mHist1D[0][i]->Fill(n);
00179 for(Int_t j = 0;j<n;j++)
00180 {
00181 StEmcCluster *c =clusters[j];
00182 if(c)
00183 {
00184 mHist1D[1][i]->Fill(c->nHits());
00185 mHist1D[2][i]->Fill(c->energy());
00186 mHist1D[3][i]->Fill(c->sigmaEta());
00187 mHist1D[4][i]->Fill(c->sigmaPhi());
00188 mHist2D[0][i]->Fill(c->eta(),c->phi());
00189 mHist2D[2][i]->Fill(c->eta(),c->phi());
00190 }
00191 }
00192 }
00193 }
00194 if(mDraw)
00195 {
00196 TString CN = "Display_";
00197 CN+=detname[i];
00198 if(!gROOT->FindObject(CN.Data()))
00199 mCanvas[i] = 0;
00200 if(!mCanvas[i])
00201 mCanvas[i] = new TCanvas(CN.Data(),CN.Data(),400,600);
00202 mCanvas[i]->cd();
00203 mHist2D[1][i]->Draw("colz");
00204 mHist2D[1][i]->GetXaxis()->UnZoom();
00205 mHist2D[1][i]->GetYaxis()->UnZoom();
00206
00207 mHist2D[2][i]->SetMarkerStyle(24);
00208 mHist2D[2][i]->SetMarkerColor(2);
00209 mHist2D[2][i]->Draw("sameP");
00210 }
00211 }
00212
00213 return StMaker::Make();
00214 }
00215
00216 Int_t StClusterDisplay::Finish()
00217 {
00218 SaveHist();
00219 return kStOk;
00220 }
00221
00222 Int_t StClusterDisplay::SaveHist()
00223 {
00224 return kStOk;
00225 }
