StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
RunGammaFitterDemo.C
1 
7 void RunGammaFitterDemo(const char* filelist = "/star/u/sdhamija/StGammaMaker/StRoot/StGammaMaker/macros/photon_9_11.list")
8 {
9  // Load shared libraries
10  gROOT->Macro("$STAR/StRoot/StGammaMaker/macros/loadGammaLibs.C");
11 
12  // Create TChain from filelist
13  TChain* chain = new TChain("gammas");
14  ifstream in(filelist);
15  string line;
16  while (getline(in, line)) chain->AddFile(line.c_str());
17  in.close();
18 
19  // Create gamma event buffer
20  StGammaEvent* event = new StGammaEvent;
21  chain->SetBranchAddress("GammaEvent", &event);
22 
23  // Create gamma fitter histograms
24  TH2* hResidualYield = new TH2F("hResidualYield", "#gamma events from pp#rightarrow#gamma+jet+X;Fit Residual Sum (R_{0,u}+R_{0,v}) [MeV];Fitted Peak Integral [MeV]", 40, 0, 200, 40, 0, 600);
25  TH1* hVertexZ = new TH1F("hVertexZ", ";z_{vertex} [cm]", 100, -200, 200);
26  TH2* hVertexXY = new TH2F("hVertexXY", ";x_{vertex} [cm];y_{vertex} [cm]", 40, -0.2, 0.2, 40, -0.45, -0.2);
27  TH1* hClusterEnergy = new TH1F("hClusterEnergy", ";Cluster E [GeV]", 100, 0, 60);
28  TH1* hClusterPt = new TH1F("hClusterPt", ";Cluster p_{T} [GeV]", 100, 0, 20);
29  TH2* hClusterXY = new TH2F("hClusterXY", ";Cluster x [cm];Cluster y [cm]", 40, -250, 250, 40, -250, 250);
30  TH2* hClusterEtaPhi = new TH2F("hClusterEtaPhi", ";Cluster #eta;Cluster #phi [radians]", 40, 1, 2, 40, -TMath::Pi(), TMath::Pi());
31  TH2* hSmdPointXY = new TH2F("hSmdPointXY", ";SMD x [cm]; SMD y [cm]", 40, -250, 250, 40, -250, 250);
32  TH1* hYieldU = new TH1F("hYieldU", ";Fitted Peak Integral [MeV]", 100, 0, 400);
33  TH1* hYieldV = new TH1F("hYieldV", ";Fitted Peak Integral [MeV]", 100, 0, 400);
34  TH1* hResidualU = new TH1F("hResidualU", ";Fit Residual [MeV]", 100, -40, 120);
35  TH1* hResidualV = new TH1F("hResidualV", ";Fit Residual [MeV]", 100, -40, 120);
36  TH1* hChi2PerNdfU = new TH1F("hChi2PerNdfU", ";#chi^{2}/ndf", 100, 0, 10);
37  TH1* hChi2PerNdfV = new TH1F("hChi2PerNdfV", ";#chi^{2}/ndf", 100, 0, 10);
38  TH1* hMeanU = new TH1F("hMeanU", ";Mean #mu", 100, 0, 288);
39  TH1* hMeanV = new TH1F("hMeanV", ";Mean #mu", 100, 0, 288);
40  TH1* hSigmaU = new TH1F("hSigmaU", ";RMS #sigma", 100, 0, 10);
41  TH1* hSigmaV = new TH1F("hSigmaV", ";RMS #sigma", 100, 0, 10);
42  TH1* hNhitsU = new TH1F("hNhitsU", ";Number of SMD hits", 50, 0, 50);
43  TH1* hNhitsV = new TH1F("hNhitsV", ";Number of SMD hits", 50, 0, 50);
44 
45  // Function for gamma/pi0 separation
46  TF1* fResidualCut = new TF1("fResidualCut", "pol2", 0, 200);
47  fResidualCut->SetParameters(100, 0, 0.05);
48 
49  // Get number of entries in gamma tree
50  int nevents = chain->GetEntries();
51  cout << "nevents = " << nevents << endl;
52 
53  // Event loop
54  int ncandidates = 0;
55  for (int iEvent = 0; iEvent < nevents; ++iEvent) {
56  chain->GetEvent(iEvent);
57  if (iEvent % 1000 == 0) cout << "iEvent = " << iEvent << endl;
58  // Skip events without vertex
59  if (event->vertex() == TVector3(0,0,0)) continue;
60  // Candidate loop
61  for (int i = 0; i < event->numberOfCandidates(); ++i) {
62  StGammaCandidate* candidate = event->candidate(i);
63  assert(candidate);
64  // EEMC candidates only
65  if (candidate->detectorId() == StGammaCandidate::kEEmc) {
66  // Fit SMD u & v plane
67  StGammaFitterResult fits[2];
68  StGammaFitterResult& u = fits[0];
69  StGammaFitterResult& v = fits[1];
70  StGammaFitter::instance()->fit(candidate, fits, 0);
71  StGammaFitter::instance()->fit(candidate, fits, 1);
72  int sector, subsector, etabin;
73  EEmcGeomSimple::Instance().getTower(candidate->position(), sector, subsector, etabin);
74  TVector3& position = EEmcSmdGeom::instance()->getIntersection(sector, u.mean, v.mean);
75  if (position.z() != -999) {
76  float residual = u.residual+v.residual;
77  float yield = u.yield+v.yield;
78  if (residual > 0 && yield > 0) {
79  ++ncandidates;
80  // Fill histograms
81  hResidualYield->Fill(residual, yield);
82  hVertexZ->Fill(event->vertex().z());
83  hVertexXY->Fill(event->vertex().x(), event->vertex().y());
84  hClusterEnergy->Fill(candidate->energy());
85  hClusterPt->Fill(candidate->momentum().Pt());
86  hClusterXY->Fill(candidate->position().x(), candidate->position().y());
87  hClusterEtaPhi->Fill(candidate->position().Eta(), candidate->position().Phi());
88  hSmdPointXY->Fill(position.x(), position.y());
89  hYieldU->Fill(u.yield);
90  hYieldV->Fill(v.yield);
91  hResidualU->Fill(u.residual);
92  hResidualV->Fill(v.residual);
93  if (u.ndf > 0) hChi2PerNdfU->Fill(u.chiSquare / u.ndf);
94  if (v.ndf > 0) hChi2PerNdfV->Fill(v.chiSquare / v.ndf);
95  hMeanU->Fill(u.mean);
96  hMeanV->Fill(v.mean);
97  hSigmaU->Fill(u.rms);
98  hSigmaV->Fill(v.rms);
99  hNhitsU->Fill(candidate->numberOfSmdu());
100  hNhitsV->Fill(candidate->numberOfSmdv());
101  }
102  }
103  }
104  } // End candidate loop
105  } // End event loop
106 
107  cout << "ncandidates = " << ncandidates << endl;
108 
109  // Create canvas
110  TCanvas* c1 = new TCanvas("c1", "c1", 1200, 800);
111  c1->Divide(3, 2);
112  gStyle->SetPalette(1);
113  gStyle->SetOptStat(11);
114 
115  // Draw histograms
116  c1->cd(1);
117  hResidualYield->Draw("colz");
118  fResidualCut->SetLineColor(kRed);
119  fResidualCut->Draw("same");
120  gPad->Print("hResidualYield.png");
121 
122  c1->cd(2);
123  hVertexZ->Draw();
124  hVertexZ->Fit("gaus");
125  hVertexZ->GetFunction("gaus")->SetLineColor(kRed);
126  gPad->Print("hVertexZ.png");
127 
128  c1->cd(3);
129  hVertexXY->Draw("colz");
130  gPad->Print("hVertexXY.png");
131 
132  c1->cd(4);
133  hClusterEnergy->Draw();
134  gPad->Print("hClusterEnergy.png");
135 
136  c1->cd(5);
137  hClusterPt->Draw();
138  gPad->Print("hClusterPt.png");
139 
140  c1->cd(6);
141  hClusterXY->Draw("colz");
142  gPad->Print("hClusterXY.png");
143 
144  c1->cd(1);
145  hClusterEtaPhi->Draw("colz");
146  gPad->Print("hClusterEtaPhi.png");
147 
148  c1->cd(2);
149  hSmdPointXY->Draw("colz");
150  gPad->Print("hSmdPointXY.png");
151 
152  c1->cd(3);
153  hYieldU->SetLineColor(kRed);
154  hYieldV->SetLineColor(kBlue);
155  hYieldU->SetLineWidth(2);
156  hYieldV->SetLineWidth(2);
157  hYieldU->Draw();
158  hYieldV->Draw("same");
159  TLegend* leg1 = new TLegend(0.60, 0.60, 0.85, 0.80);
160  leg1->AddEntry(hYieldU, "SMD-u", "L");
161  leg1->AddEntry(hYieldV, "SMD-v", "L");
162  leg1->Draw();
163  gPad->Print("hYield.png");
164 
165  c1->cd(4);
166  hResidualU->SetLineColor(kRed);
167  hResidualV->SetLineColor(kBlue);
168  hResidualU->SetLineWidth(2);
169  hResidualV->SetLineWidth(2);
170  hResidualU->Draw();
171  hResidualV->Draw("same");
172  TLegend* leg2 = new TLegend(0.60, 0.60, 0.85, 0.80);
173  leg2->AddEntry(hResidualU, "SMD-u", "L");
174  leg2->AddEntry(hResidualV, "SMD-v", "L");
175  leg2->Draw();
176  gPad->Print("hResidual.png");
177 
178  c1->cd(5);
179  hChi2PerNdfU->SetLineColor(kRed);
180  hChi2PerNdfV->SetLineColor(kBlue);
181  hChi2PerNdfU->SetLineWidth(2);
182  hChi2PerNdfV->SetLineWidth(2);
183  hChi2PerNdfU->Draw();
184  hChi2PerNdfV->Draw("same");
185  TLegend* leg3 = new TLegend(0.60, 0.60, 0.85, 0.80);
186  leg3->AddEntry(hChi2PerNdfU, "SMD-u", "L");
187  leg3->AddEntry(hChi2PerNdfV, "SMD-v", "L");
188  leg3->Draw();
189  gPad->Print("hChi2PerNdf.png");
190 
191  c1->cd(6);
192  hMeanU->SetLineColor(kRed);
193  hMeanV->SetLineColor(kBlue);
194  hMeanU->SetLineWidth(2);
195  hMeanV->SetLineWidth(2);
196  hMeanU->Draw();
197  hMeanV->Draw("same");
198  TLegend* leg4 = new TLegend(0.60, 0.60, 0.85, 0.80);
199  leg4->AddEntry(hMeanU, "SMD-u", "L");
200  leg4->AddEntry(hMeanV, "SMD-v", "L");
201  leg4->Draw();
202  gPad->Print("hMean.png");
203 
204  c1->cd(1);
205  hSigmaU->SetLineColor(kRed);
206  hSigmaV->SetLineColor(kBlue);
207  hSigmaU->SetLineWidth(2);
208  hSigmaV->SetLineWidth(2);
209  hSigmaU->Draw();
210  hSigmaV->Draw("same");
211  TLegend* leg5 = new TLegend(0.60, 0.60, 0.85, 0.80);
212  leg5->AddEntry(hSigmaU, "SMD-u", "L");
213  leg5->AddEntry(hSigmaV, "SMD-v", "L");
214  leg5->Draw();
215  gPad->Print("hSigma.png");
216 
217  c1->cd(2);
218  hNhitsU->SetLineColor(kRed);
219  hNhitsV->SetLineColor(kBlue);
220  hNhitsU->SetLineWidth(2);
221  hNhitsV->SetLineWidth(2);
222  hNhitsU->Draw();
223  hNhitsV->Draw("same");
224  TLegend* leg6 = new TLegend(0.60, 0.60, 0.85, 0.80);
225  leg6->AddEntry(hNhitsU, "SMD-u", "L");
226  leg6->AddEntry(hNhitsV, "SMD-v", "L");
227  leg6->Draw();
228  gPad->Print("hNhits.png");
229 }
EEmcGeomSimple::Instance
static EEmcGeomSimple & Instance()
returns a reference to a static instance of EEmcGeomSimple
Definition: EEmcGeomSimple.h:178
StGammaFitter::instance
static StGammaFitter * instance()
Access to single instance of this singleton class.
Definition: StGammaFitter.cxx:67
StGammaFitter::fit
int fit(StGammaCandidate *candidate, StGammaFitterResult *fits, Int_t plane=0)
Fit transverse SMD profile to predetermined peak in u- and v-plane.
Definition: StGammaFitter.cxx:73
EEmcSmdGeom::getIntersection
TVector3 getIntersection(Int_t iSec, Int_t iUStrip, Int_t iVStrip, const TVector3 &vertex) const
Definition: EEmcSmdGeom.cxx:628
EEmcGeomSimple::getTower
bool getTower(const TVector3 &r, int &sec, int &sub, int &etabin, Float_t &dphi, Float_t &deta) const
Definition: EEmcGeomSimple.cxx:159
StGammaEvent::candidate
StGammaCandidate * candidate(Int_t i) const
Return ith strip.
Definition: StGammaEvent.h:71
Generated by   doxygen 1.8.5