StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StMcQaMaker.cxx
1 /*************************************************
2  *
3  * $Id: StMcQaMaker.cxx,v 1.4 2010/06/01 20:46:47 perev Exp $
4  * Author: Manuel Calderon de la Barca
5  * Make standard Histograms for
6  * StMcEvent
7  * $Log: StMcQaMaker.cxx,v $
8  * Revision 1.4 2010/06/01 20:46:47 perev
9  * const added
10  *
11  * Revision 1.3 2010/01/28 18:13:18 perev
12  * WarningOff
13  *
14  * Revision 1.2 2007/10/17 19:32:06 fisyak
15  * The pixel story, remove Igt and Fst
16  *
17  * Revision 1.1 2007/03/21 16:48:49 fisyak
18  * maker for side by side comparision of GEANT and VMC simulations
19  *
20  *
21  *************************************************/
22 //#include <assert.h>
23 #include <Stiostream.h>
24 #include <stdlib.h>
25 //#include <string>
26 //#include <vector>
27 //#include <cmath>
28 
29 #include "TH1.h"
30 #include "TH2.h"
31 #include "TFile.h"
32 
33 #include "StMcQaMaker.h"
34 //#include "PhysicalConstants.h"
35 //#include "SystemOfUnits.h"
36 
37 #include "StThreeVectorF.hh"
38 #include "StEnumerations.h"
39 #include "StMcEventTypes.hh"
40 
41 #include "StMcEvent.hh"
42 
43 #include "StEmcUtil/geometry/StEmcGeom.h"
44 #include "TMath.h"
45 ClassImp(StMcQaMaker)
46 
47 //_________________________________________________
48 Int_t StMcQaMaker::Init() {
49  QAPlots(0);
50  return StMaker::Init();
51 }
52 //_________________________________________________
53 Int_t StMcQaMaker::Make() {
54  // StMcEvent
55  StMcEvent* mEvent = (StMcEvent*) GetDataSet("StMcEvent");
56  if (mEvent) QAPlots(mEvent);
57  return kStOk;
58 }
59 
60 //--------------------------------------------------------------------------------
61 void StMcQaMaker::QAPlots(StMcEvent *mEvent) {
62  // Histograms of StMcEvent
63  // quantities.
64  struct TrackDetector_t {
65  StDetectorId Id;
66  const Char_t *Name;
67  Int_t NHMax; // max. no. of hit per track per detector
68  TH1F *NHits[2];//
69  Int_t NZ;
70  Double_t Zmax;
71  TH1F *Zhit[2];
72  Int_t NX, NY;
73  Double_t Xmax, Ymax;
74  TH2F *XYhit[2];
75  TH2F *TdEdxP[2];
76  };
77  struct CalorimeterDetector_t {
78  StDetectorId Id;
79  const Char_t *Name;
80  Int_t NHMax; // max. no. of hit per track per detector
81  TH1F *NHits[2];//
82  Int_t Neta, Nphi;
83  Double_t EtaMin,EtaMax;
84  TH2F *EtaPhiHit[2];
85  TH2F *dEEtaPhiHit[2];
86  };
87  static TrackDetector_t TrackDetectors[] = {
88  // Id Name NHMax NZ NX NY XY dEdx
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}},
92  // {kFtpcWestId ,"FtpcWest", 10, {0, 0}, 60, 60., {0, 0}, 100, 100, 20., 20., {0, 0}, {0, 0}},
93  // {kFtpcEastId ,"FtpcEast", 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}},
98  // {kZdcWestId ,"ZdcWest", 10, {0, 0}, 250, 250., {0, 0}, 250, 250, 250., 250., {0, 0}, {0, 0}},
99  // {kZdcEastId ,"ZdcEast", 10, {0, 0}, 250, 250., {0, 0}, 250, 250, 250., 250., {0, 0}, {0, 0}},
100  {kZdcEastId ,"Zdc", 10, {0, 0}, 250, 250., {0, 0}, 250, 250, 250., 250., {0, 0}, {0, 0}},
101  // {kMwpcWestId ,"MwpcWest", 10, {0, 0}, 250, 250., {0, 0}, 250, 250, 250., 250., {0, 0}, {0, 0}},
102  // {kMwpcEastId ,"MwpcEast", 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}}
109  };
110  static const Int_t NTdetectors = sizeof(TrackDetectors)/sizeof(TrackDetector_t);
111  static CalorimeterDetector_t CalorimeterDetectors[] = {
112  // Id Name NHMax Neta Nphi EtaMin Max EtaPhi dEEtaPhi
113  {kBarrelEmcTowerId ,"Bemc", 10, {0, 0}, 40, 480, -1., 1., {0, 0}, {0, 0}}, // (D_eta,D_phi) = (0.05,0.05)
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}}
121  };
122  static const Int_t NCdetectors = sizeof(CalorimeterDetectors)/sizeof(CalorimeterDetector_t);
123  // Detector Information
124  // per-Event [k] = 0 for primaries, 1 for all
125  static TH1F* mVertexZ[2] = {0,0};
126  static TH2F* mVertexXY[2] = {0,0};
127  // Track Information
128 
129  // eventwise particle multiplicity
130  static TH1F* mTracks[2] = {0, 0};
131  static TH1F* mGeantIds[2] = {0, 0};
132 
133  // trackwise distributions
134  static TH1F* mPt[2] = {0, 0};
135  static TH1F* mEta[2] = {0, 0};
136  if (! mEvent) {// Book
137  Int_t Begin = gDirectory->GetList()->GetSize();
138  TString Names[2] = {"Primary","All"};
139  TString Name("");
140  TString Title("");
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);
154  // Track detectors
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);
171  }
172  // Calorimeters
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);
187  }
188  }
189  // Move them into hist TDataSet
190  Int_t Last = gDirectory->GetList()->GetSize() -1;
191  for (Int_t i = Last; i >= Begin; i--) AddHist((TH1*) gDirectory->GetList()->At(i));
192  return;
193  }
194  // Fill
195  //
196  // Vertex
197  //
198  for (Int_t k = 0; k < 2; k++) {
199  StMcVertex *vertex = 0;
200  StPtrVecMcTrack *tracks = 0;
201  Int_t nVertices = 0;
202  if (k == 0) {
203  vertex = mEvent->primaryVertex();
204  nVertices = 1;
205  tracks = &vertex->daughters();
206  } else {
207  nVertices = mEvent->vertices().size();
208  tracks = &mEvent->tracks();
209  }
210  for (int i = 0; i < nVertices; i++) {
211  if (k) vertex = (mEvent->vertices()[i]);
212  if (vertex) {
213  mVertexZ[k]->Fill(vertex->position().z());
214  mVertexXY[k]->Fill(vertex->position().x(),vertex->position().y());
215  if (tracks) {
216  StPtrVecMcTrack &Tracks = *tracks;
217  mTracks[k]->Fill(Tracks.size());
218  for (size_t iTrk=0; iTrk<Tracks.size(); ++iTrk) {
219  StMcTrack* trk = Tracks[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);
224  // hits
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++) {
232  StMcHit *hit = Hits[h];
233  if (! hit) continue;
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;
239  dEdx /= hit->dS();
240  Double_t pmom = hit->localMomentum().mag();
241  TrackDetectors[det].TdEdxP[k]->Fill(TMath::Log10(pmom),TMath::Log10(dEdx));
242  }
243  }
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++) {
251  StMcCalorimeterHit *hit = Hits[h];
252  if (! hit) continue;
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());
255  }
256  }
257  }
258  }
259  }
260  }
261  }
262 }
StMcQaMaker::Make
virtual Int_t Make()
Definition: StMcQaMaker.cxx:53
TDataSet::Last
virtual TDataSet * Last() const
Return the last object in the list. Returns 0 when list is empty.
Definition: TDataSet.cxx:437
StMcTrack
Monte Carlo Track class All information on a simulated track is stored in this class: kinematics...
Definition: StMcTrack.hh:144
StMcQaMaker::QAPlots
virtual void QAPlots(StMcEvent *ev=0)
Definition: StMcQaMaker.cxx:61
StMcEvent
Event data structure to hold all information from a Monte Carlo simulation. This class is the interfa...
Definition: StMcEvent.hh:169
kStOk
Definition: Stypes.h:41
Generated by   doxygen 1.8.5