33 #include "StMcQaMaker.h"
37 #include "StThreeVectorF.hh"
39 #include "StMcEventTypes.hh"
41 #include "StMcEvent.hh"
43 #include "StEmcUtil/geometry/StEmcGeom.h"
50 return StMaker::Init();
64 struct TrackDetector_t {
77 struct CalorimeterDetector_t {
83 Double_t EtaMin,EtaMax;
87 static TrackDetector_t TrackDetectors[] = {
89 {kTpcId ,
"Tpc", 50, {0, 0}, 210, 210., {0, 0}, 100, 100, 200., 200., {0, 0}, {0, 0}},
90 {kSvtId ,
"Svt", 10, {0, 0}, 30, 30., {0, 0}, 100, 100, 20., 20., {0, 0}, {0, 0}},
91 {kRichId ,
"Rich", 10, {0, 0}, 60, 60., {0, 0}, 100, 100, 20., 20., {0, 0}, {0, 0}},
94 {kFtpcEastId ,
"Ftpc", 10, {0, 0}, 150, 300., {0, 0}, 100, 100, 40., 40., {0, 0}, {0, 0}},
95 {kTofId ,
"Tof", 10, {0, 0}, 200, 200., {0, 0}, 100, 100, 200., 200., {0, 0}, {0, 0}},
96 {kCtbId ,
"Ctb", 10, {0, 0}, 100, 300., {0, 0}, 100, 100, 300., 300., {0, 0}, {0, 0}},
97 {kSsdId ,
"Ssd", 10, {0, 0}, 60, 60., {0, 0}, 100, 100, 40., 40., {0, 0}, {0, 0}},
100 {kZdcEastId ,
"Zdc", 10, {0, 0}, 250, 250., {0, 0}, 250, 250, 250., 250., {0, 0}, {0, 0}},
103 {kMwpcEastId ,
"Mwpc", 10, {0, 0}, 250, 250., {0, 0}, 250, 250, 250., 250., {0, 0}, {0, 0}},
104 {kPhmdCpvId ,
"PhmdCpv", 10, {0, 0}, 250, 250., {0, 0}, 250, 250, 250., 250., {0, 0}, {0, 0}},
105 {kPhmdId ,
"Phmd", 10, {0, 0}, 250, 250., {0, 0}, 250, 250, 250., 250., {0, 0}, {0, 0}},
106 {kPxlId ,
"Pxl", 10, {0, 0}, 60, 60., {0, 0}, 100, 100, 20., 20., {0, 0}, {0, 0}},
107 {kIstId ,
"Ist", 10, {0, 0}, 60, 60., {0, 0}, 100, 100, 20., 20., {0, 0}, {0, 0}},
108 {kFgtId ,
"Fgt", 10, {0, 0}, 250, 250., {0, 0}, 250, 250, 250., 250., {0, 0}, {0, 0}}
110 static const Int_t NTdetectors =
sizeof(TrackDetectors)/
sizeof(TrackDetector_t);
111 static CalorimeterDetector_t CalorimeterDetectors[] = {
113 {kBarrelEmcTowerId ,
"Bemc", 10, {0, 0}, 40, 480, -1., 1., {0, 0}, {0, 0}},
114 {kBarrelEmcPreShowerId ,
"Bprs", 10, {0, 0}, 240, 240, -1., 1., {0, 0}, {0, 0}},
115 {kBarrelSmdEtaStripId ,
"Bsmde", 10, {0, 0}, 240, 480, -1., 1., {0, 0}, {0, 0}},
116 {kBarrelSmdPhiStripId ,
"Bsmdp", 10, {0, 0}, 40, 480, -1., 1., {0, 0}, {0, 0}},
117 {kEndcapEmcTowerId ,
"Eemc", 10, {0, 0}, 120, 240, 1., 3., {0, 0}, {0, 0}},
118 {kEndcapEmcPreShowerId ,
"Eprs", 10, {0, 0}, 120, 240, 1., 3., {0, 0}, {0, 0}},
119 {kEndcapSmdUStripId ,
"Esmdu", 10, {0, 0}, 120, 240, 1., 3., {0, 0}, {0, 0}},
120 {kEndcapSmdVStripId ,
"Esmdv", 10, {0, 0}, 120, 240, 1., 3., {0, 0}, {0, 0}}
122 static const Int_t NCdetectors =
sizeof(CalorimeterDetectors)/
sizeof(CalorimeterDetector_t);
125 static TH1F* mVertexZ[2] = {0,0};
126 static TH2F* mVertexXY[2] = {0,0};
130 static TH1F* mTracks[2] = {0, 0};
131 static TH1F* mGeantIds[2] = {0, 0};
134 static TH1F* mPt[2] = {0, 0};
135 static TH1F* mEta[2] = {0, 0};
137 Int_t Begin = gDirectory->GetList()->GetSize();
138 TString Names[2] = {
"Primary",
"All"};
141 for (Int_t k = 0; k < 2; k++) {
142 Name =
"McQAVertexZ"; Name += Names[k];
143 mVertexZ[k] =
new TH1F(Name,
"Vertex Z Position",300,-300,300);
144 Name =
"McQAVertexXY"; Name += Names[k];
145 mVertexXY[k] =
new TH2F(Name,
"Vertex XY Position",200,-300,300,200,-300,300);
146 Name =
"McQATracks"; Name += Names[k];
147 mTracks[k] =
new TH1F(Name,
"All Tracks",1000,0,10000);
148 Name =
"McGeantIds"; Name += Names[k];
149 mGeantIds[k] =
new TH1F(Name,
"GeantId for all tracks", 50, 0, 50);
150 Name =
"McQApT"; Name += Names[k];
151 mPt[k] =
new TH1F(Name,
"pT",200,0,10);
152 Name =
"McQAEta"; Name += Names[k];
153 mEta[k] =
new TH1F(Name,
"eta",200,-5,5);
155 for (Int_t det = 0; det < NTdetectors; det++) {
156 Name =
"McQAN"; Name += TrackDetectors[det].Name; Name += Names[k];
157 Title =
"No. of Hits in "; Title += TrackDetectors[det].Name; Title +=
" for "; Title += Names[k]; Title +=
" tracks";
158 TrackDetectors[det].NHits[k] =
new TH1F(Name,Title,TrackDetectors[det].NHMax,0,TrackDetectors[det].NHMax);
159 Name =
"McQAZ"; Name += TrackDetectors[det].Name; Name += Names[k];
160 Title =
"Z of Hit in "; Title += TrackDetectors[det].Name; Title +=
" for "; Title += Names[k]; Title +=
" tracks";
161 TrackDetectors[det].Zhit[k] =
new TH1F(Name,Title,TrackDetectors[det].NZ,-TrackDetectors[det].Zmax,TrackDetectors[det].Zmax);
162 Name =
"McQAXY"; Name += TrackDetectors[det].Name; Name += Names[k];
163 Title =
"XY of Hit in "; Title += TrackDetectors[det].Name; Title +=
" for "; Title += Names[k]; Title +=
" tracks";
164 TrackDetectors[det].XYhit[k] =
new TH2F(Name,Title,
165 TrackDetectors[det].NX,-TrackDetectors[det].Xmax,TrackDetectors[det].Xmax,
166 TrackDetectors[det].NY,-TrackDetectors[det].Ymax,TrackDetectors[det].Ymax);
167 Name =
"McQAdEdx"; Name += TrackDetectors[det].Name; Name += Names[k];
168 Title =
"log10(dE/dx) (keV/cm)) versus log10(p(GeV/c) in ";
169 Title += TrackDetectors[det].Name; Title +=
" for "; Title += Names[k]; Title +=
" tracks";
170 TrackDetectors[det].TdEdxP[k] =
new TH2F(Name,Title,200,-2.,2., 1600,-1.,7.0);
173 for (Int_t det = 0; det < NCdetectors; det++) {
174 Name =
"McQAN"; Name += CalorimeterDetectors[det].Name; Name += Names[k];
175 Title =
"No. of Hits in "; Title += CalorimeterDetectors[det].Name; Title +=
" for "; Title += Names[k]; Title +=
" tracks";
176 CalorimeterDetectors[det].NHits[k] =
new TH1F(Name,Title,CalorimeterDetectors[det].NHMax,0,CalorimeterDetectors[det].NHMax);
177 Name =
"McQAEP"; Name += CalorimeterDetectors[det].Name; Name += Names[k];
178 Title =
"Eta/Phi of Hit in "; Title += CalorimeterDetectors[det].Name; Title +=
" for "; Title += Names[k]; Title +=
" tracks";
179 CalorimeterDetectors[det].EtaPhiHit[k] =
new TH2F(Name,Title,
180 CalorimeterDetectors[det].Neta,0,CalorimeterDetectors[det].Neta,
181 CalorimeterDetectors[det].Nphi,0,CalorimeterDetectors[det].Nphi);
182 Name =
"McQAdEEP"; Name += CalorimeterDetectors[det].Name; Name += Names[k];
183 Title =
"Sum of dE"; Title += CalorimeterDetectors[det].Name; Title +=
" for "; Title += Names[k]; Title +=
" tracks";
184 CalorimeterDetectors[det].dEEtaPhiHit[k] =
new TH2F(Name,Title,
185 CalorimeterDetectors[det].Neta,0,CalorimeterDetectors[det].Neta,
186 CalorimeterDetectors[det].Nphi,0,CalorimeterDetectors[det].Nphi);
190 Int_t
Last = gDirectory->GetList()->GetSize() -1;
191 for (Int_t i = Last; i >= Begin; i--) AddHist((TH1*) gDirectory->GetList()->At(i));
198 for (Int_t k = 0; k < 2; k++) {
200 StPtrVecMcTrack *tracks = 0;
203 vertex = mEvent->primaryVertex();
205 tracks = &vertex->daughters();
207 nVertices = mEvent->vertices().size();
208 tracks = &mEvent->tracks();
210 for (
int i = 0; i < nVertices; i++) {
211 if (k) vertex = (mEvent->vertices()[i]);
213 mVertexZ[k]->Fill(vertex->position().z());
214 mVertexXY[k]->Fill(vertex->position().x(),vertex->position().y());
216 StPtrVecMcTrack &Tracks = *tracks;
217 mTracks[k]->Fill(Tracks.size());
218 for (
size_t iTrk=0; iTrk<Tracks.size(); ++iTrk) {
220 mPt[k]->Fill(trk->momentum().perp());
221 mEta[k]->Fill(trk->momentum().pseudoRapidity());
222 Int_t gId = trk->geantId();
223 mGeantIds[k]->Fill(gId);
225 for (Int_t det = 0; det < NTdetectors; det++) {
226 const StPtrVecMcHit *hits = trk->Hits(TrackDetectors[det].Id);
227 if (! hits)
continue;
228 const StPtrVecMcHit &Hits = *hits;
229 Int_t Nhits = Hits.size();
230 TrackDetectors[det].NHits[k]->Fill(Nhits);
231 for (Int_t h = 0; h < Nhits; h++) {
234 TrackDetectors[det].Zhit[k]->Fill(hit->position().z());
235 TrackDetectors[det].XYhit[k]->Fill(hit->position().x(),hit->position().y());
236 Double_t dEdx = 1e6*TMath::Abs(hit->dE());
237 if (dEdx < 1.e-10)
continue;
238 if (hit->dS() < 1.e-10)
continue;
240 Double_t pmom = hit->localMomentum().mag();
241 TrackDetectors[det].TdEdxP[k]->Fill(TMath::Log10(pmom),TMath::Log10(dEdx));
244 for (Int_t det = 0; det < NCdetectors; det++) {
245 const StPtrVecMcCalorimeterHit *hits = trk->CalorimeterHits(CalorimeterDetectors[det].Id);
246 if (! hits)
continue;
247 const StPtrVecMcCalorimeterHit &Hits = *hits;
248 Int_t Nhits = Hits.size();
249 CalorimeterDetectors[det].NHits[k]->Fill(Nhits);
250 for (Int_t h = 0; h < Nhits; h++) {
253 CalorimeterDetectors[det].EtaPhiHit[k]->Fill(hit->eta(),2*hit->module()+hit->sub());
254 CalorimeterDetectors[det].dEEtaPhiHit[k]->Fill(hit->eta(),2*hit->module()+hit->sub(),hit->dE());
virtual TDataSet * Last() const
Return the last object in the list. Returns 0 when list is empty.
Monte Carlo Track class All information on a simulated track is stored in this class: kinematics...
virtual void QAPlots(StMcEvent *ev=0)
Event data structure to hold all information from a Monte Carlo simulation. This class is the interfa...