Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
OccupancySim.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file OccupancySim.C
1 
2 
3 #include "SaveCanvas.C"
4 #include "sPhenixStyle.C"
5 
6 #include <TCanvas.h>
7 #include <TDatabasePDG.h>
8 #include <TFile.h>
9 #include <TH1D.h>
10 #include <TH1F.h>
11 #include <TH2F.h>
12 #include <TH3F.h>
13 #include <TLegend.h>
14 #include <TLine.h>
15 #include <TLorentzVector.h>
16 #include <TRandom3.h>
17 #include <TString.h>
18 #include <TTree.h>
19 #include <TVirtualFitter.h>
20 
21 #include <assert.h>
22 
23 #include <cmath>
24 #include <cstddef>
25 #include <iostream>
26 
27 using namespace std;
28 
29 TFile *_file0 = NULL;
30 TFile *_file_trigger = NULL;
31 TString description;
32 
33 typedef vector<vector<TH1 *> > chipMultiplicitySet_vec;
34 
35 chipMultiplicitySet_vec MakeChipMultiplicitySet(TFile *file)
36 {
37  assert(file);
38 
39  TH3 *LayerChipMultiplicity = (TH3 *) file->GetObjectChecked("LayerChipMultiplicity", "TH3");
40 
41  TCanvas *c1 = new TCanvas("MakeChipMultiplicitySet", "MakeChipMultiplicitySet", 1800, 960);
42  c1->Divide(9, 3, 0, 0);
43  int idx = 1;
44  TPad *p;
45 
46  chipMultiplicitySet_vec chipMultiplicitySet;
47 
48  for (int ilayer = 1; ilayer <= LayerChipMultiplicity->GetNbinsX(); ++ilayer)
49  {
50  vector<TH1 *> chipMultiplicityLayer;
51 
52  for (int ichip = 1; ichip <= LayerChipMultiplicity->GetNbinsY(); ++ichip)
53  {
54  p = (TPad *) c1->cd(idx++);
55  c1->Update();
56  p->SetLogy();
57 
58  LayerChipMultiplicity->GetXaxis()->SetRange(ilayer, ilayer);
59  LayerChipMultiplicity->GetYaxis()->SetRange(ichip, ichip);
60 
61  TString histname(Form("ChipMultiplicity_Layer%dChip%d", ilayer - 1, ichip - 1));
62  TH1 *h = LayerChipMultiplicity->Project3D("z");
63  h->SetTitle(";Chip Multiplicity [pixel];Count");
64 
65  h->SetName(histname);
66 
67  h->Draw();
68 
69  h->SetDirectory(NULL);
70  h->ComputeIntegral(true);
71 
72  chipMultiplicityLayer.push_back(h);
73  }
74 
75  chipMultiplicitySet.push_back(chipMultiplicityLayer);
76  }
77 
78  SaveCanvas(c1, TString(file->GetName()) + TString(c1->GetName()), false);
79 
80  c1 = new TCanvas("MakeChipMultiplicitySetChip4", "MakeChipMultiplicitySetChip4", 1200, 960);
81  gPad->SetLogy();
82  TH1 *LayerMultiplicityLayer0 = chipMultiplicitySet[0][4];
83  TH1 *LayerMultiplicityLayer1 = chipMultiplicitySet[1][4];
84  TH1 *LayerMultiplicityLayer2 = chipMultiplicitySet[2][4];
85 
86  LayerMultiplicityLayer0->SetLineColor(kRed + 2);
87  LayerMultiplicityLayer1->SetLineColor(kBlue + 2);
88  LayerMultiplicityLayer2->SetLineColor(kGreen + 2);
89 
90  LayerMultiplicityLayer0->Draw();
91  LayerMultiplicityLayer1->Draw("same");
92  LayerMultiplicityLayer2->Draw("same");
93 
94  LayerMultiplicityLayer0->SetTitle(";Chip multiplicity [Pixel];Count");
95 
96  TLegend *leg = new TLegend(.5, .5, .93, .93);
97  leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation", "");
98  leg->AddEntry("", description, "");
99  leg->AddEntry(LayerMultiplicityLayer0, Form("MVTX Layer0 Chip4, <hit> = %.1f", LayerMultiplicityLayer0->GetMean()), "l");
100  leg->AddEntry(LayerMultiplicityLayer1, Form("MVTX Layer1 Chip4, <hit> = %.1f", LayerMultiplicityLayer1->GetMean()), "l");
101  leg->AddEntry(LayerMultiplicityLayer2, Form("MVTX Layer2 Chip4, <hit> = %.1f", LayerMultiplicityLayer2->GetMean()), "l");
102  leg->Draw();
103 
104  SaveCanvas(c1, TString(file->GetName()) + TString(c1->GetName()), false);
105 
106  return chipMultiplicitySet;
107 }
108 
109 void Check(TFile *file)
110 {
111  assert(file);
112 
113  TCanvas *c1 = new TCanvas("Check", "Check", 1800, 960);
114  c1->Divide(2, 2);
115  int idx = 1;
116  TPad *p;
117 
118  p = (TPad *) c1->cd(idx++);
119  c1->Update();
120  // p->SetLogz();
121 
122  TH1 *hNormalization = (TH1 *) file->GetObjectChecked("hNormalization", "TH1");
123  hNormalization->Draw();
124 
125  p = (TPad *) c1->cd(idx++);
126  c1->Update();
127  // p->SetLogy();
128 
129  TH1 *hNChEta = (TH1 *) file->GetObjectChecked("hNChEta", "TH1");
130  assert(hNChEta);
131  double norm = hNChEta->GetBinWidth(1) * hNormalization->GetBinContent(2);
132  hNChEta->Scale(1. / norm);
133  hNChEta->Draw();
134 
135  hNChEta->GetYaxis()->SetTitle("dN_{Ch}/d#eta");
136  // hNormalization->Draw()
137 
138  p = (TPad *) c1->cd(idx++);
139  c1->Update();
140  TH1 *hVertexZ = (TH1 *) file->GetObjectChecked("hVertexZ", "TH1");
141  hVertexZ->Draw();
142 
143  p = (TPad *) c1->cd(idx++);
144  c1->Update();
145  p->SetLogy();
146 
147  TH2 *LayerMultiplicity = (TH2 *) file->GetObjectChecked("LayerMultiplicity", "TH2");
148  TH1 *LayerMultiplicityLayer0 =
149  LayerMultiplicity->ProjectionY("LayerMultiplicityLayer0", 1, 1);
150  TH1 *LayerMultiplicityLayer1 =
151  LayerMultiplicity->ProjectionY("LayerMultiplicityLayer1", 2, 2);
152  TH1 *LayerMultiplicityLayer2 =
153  LayerMultiplicity->ProjectionY("LayerMultiplicityLayer2", 3, 3);
154 
155  LayerMultiplicityLayer0->SetLineColor(kRed + 2);
156  LayerMultiplicityLayer1->SetLineColor(kBlue + 2);
157  LayerMultiplicityLayer2->SetLineColor(kGreen + 2);
158 
159  LayerMultiplicityLayer0->Draw();
160  LayerMultiplicityLayer1->Draw("same");
161  LayerMultiplicityLayer2->Draw("same");
162 
163  LayerMultiplicityLayer0->SetTitle(";Chip multiplicity [Pixel];Count");
164 
165  TLegend *leg = new TLegend(.5, .5, .93, .93);
166  leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation", "");
167  leg->AddEntry("", description, "");
168  leg->AddEntry(LayerMultiplicityLayer0, Form("MVTX Layer0, <hit> = %.1f", LayerMultiplicityLayer0->GetMean()), "l");
169  leg->AddEntry(LayerMultiplicityLayer1, Form("MVTX Layer1, <hit> = %.1f", LayerMultiplicityLayer1->GetMean()), "l");
170  leg->AddEntry(LayerMultiplicityLayer2, Form("MVTX Layer2, <hit> = %.1f", LayerMultiplicityLayer2->GetMean()), "l");
171  leg->Draw();
172 
173  SaveCanvas(c1, TString(file->GetName()) + TString(c1->GetName()), false);
174 }
175 
176 TH1 *MakeCDF(TH1 *h)
177 {
178  assert(h);
179 
180  TH1 *hCDF = (TH1 *) h->Clone(TString("CDF") + h->GetName());
181  hCDF->SetDirectory(NULL);
182 
183  double integral = 0;
184 
185  for (int bin = h->GetNbinsX() + 1; bin >= 0; --bin)
186  {
187  integral += h->GetBinContent(bin);
188  hCDF->SetBinContent(bin, integral);
189  }
190 
191  for (int bin = h->GetNbinsX(); bin >= 1; --bin)
192  {
193  hCDF->SetBinContent(bin, hCDF->GetBinContent(bin) / integral);
194  }
195 
196  return hCDF;
197 }
198 
199 chipMultiplicitySet_vec TriggerMultiplicity(chipMultiplicitySet_vec cm_MB, chipMultiplicitySet_vec cm_Trigger, const double mu_MB, const int n_Trigger, const double mu_Noise)
200 {
201  //init
202  chipMultiplicitySet_vec cm;
203 
204  const static int MaxMultiplicity(2000);
205 
206  for (auto &layer : cm_MB)
207  {
208  vector<TH1 *> chipMultiplicityLayer;
209 
210  for (auto &chip : layer)
211  {
212  TH1 *h = new TH1D(
213  TString("Trigger") + chip->GetName(), TString("Trigger") + chip->GetName(),
214  MaxMultiplicity, -.5, MaxMultiplicity - .5);
215  // (TH1 *) chip->Clone(TString("Trigger") + chip->GetName());
216 
217  h->SetDirectory(NULL);
218  h->GetXaxis()->SetTitle(chip->GetXaxis()->GetTitle());
219  h->GetYaxis()->SetTitle(chip->GetYaxis()->GetTitle());
220 
221  assert(h);
222  chipMultiplicityLayer.push_back(h);
223  }
224  cm.push_back(chipMultiplicityLayer);
225  }
226 
227  //composition for trigger
228  TRandom3 rnd(1234);
229  assert(mu_MB >= 0);
230  assert(n_Trigger >= 0 && n_Trigger <= 1);
231  assert(mu_Noise >= 0);
232 
233  // const int NSample = 1000000;
234  const int NSample = 10000000;
235  for (int i = 0; i < NSample; ++i)
236  {
237  if (i % (NSample / 10) == 0)
238  {
239  cout << "TriggerMultiplicity - " << i << endl;
240  }
241 
242  int n_MB = rnd.Poisson(mu_MB);
243  int n_Noise = rnd.Poisson(mu_Noise);
244 
245  int ilayer = 0;
246  for (auto &layer : cm)
247  {
248  int ichip = 0;
249  for (auto &chip : layer)
250  {
251  assert(chip);
252  int n_hit = n_Noise;
253 
254  assert(ilayer < cm_Trigger.size());
255  assert(ichip < cm_Trigger[ilayer].size());
256  assert(cm_Trigger[ilayer][ichip]);
257  if (n_Trigger)
258  {
259  double rndHit = cm_Trigger[ilayer][ichip]->GetRandom();
260  if (rndHit < .5) rndHit = 0;
261 
262  // cout << "n_hit += " << cm_Trigger[ilayer][ichip]->GetRandom() << endl;
263  n_hit += rndHit;
264  }
265 
266  assert(ilayer < cm_MB.size());
267  assert(ichip < cm_MB[ilayer].size());
268  assert(cm_MB[ilayer][ichip]);
269  for (int iMB = 0; iMB < n_MB; ++iMB)
270  {
271  double rndHit = cm_MB[ilayer][ichip]->GetRandom();
272  if (rndHit < .5) rndHit = 0;
273  // cout << "n_hit +=" << cm_MB[ilayer][ichip]->GetRandom() << endl;
274  n_hit += rndHit;
275  }
276 
277  // cout << "chip->Fill(n_hit);" << n_hit << ", n_Noise = " << n_Noise << endl;
278  chip->Fill(n_hit);
279 
280  ++ichip;
281  }
282  ++ilayer;
283  }
284  }
285 
286  TString UniqueName(Form("_muMB%.0f_nTrig%d_muNoise%.0f", mu_MB * 100, n_Trigger, mu_Noise));
287 
288  //save - plot
289  TCanvas *c1 = new TCanvas("TriggerMultiplicity" + UniqueName, "TriggerMultiplicity" + UniqueName, 1800, 960);
290  c1->Divide(9, 3, 0, 0);
291  int idx = 1;
292  TPad *p;
293 
294  for (auto &layer : cm)
295  {
296  vector<TH1 *> chipMultiplicityLayer;
297 
298  for (auto &chip : layer)
299  {
300  p = (TPad *) c1->cd(idx++);
301  c1->Update();
302  p->SetLogy();
303 
304  chip->SetTitle(";Chip Multiplicity [pixel];Count");
305  chip->Draw();
306  }
307  }
308 
309  SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), false);
310 
311  c1 = new TCanvas("TriggerMultiplicityChip4" + UniqueName, "TriggerMultiplicityChip4" + UniqueName, 1900, 750);
312  c1->Divide(2, 1);
313  idx = 1;
314 
315  p = (TPad *) c1->cd(idx++);
316  c1->Update();
317  p->SetLogy();
318 
319  TH1 *LayerMultiplicityLayer0 = cm[0][4];
320  TH1 *LayerMultiplicityLayer1 = cm[1][4];
321  TH1 *LayerMultiplicityLayer2 = cm[2][4];
322 
323  LayerMultiplicityLayer0->SetLineColor(kRed + 2);
324  LayerMultiplicityLayer1->SetLineColor(kBlue + 2);
325  LayerMultiplicityLayer2->SetLineColor(kGreen + 2);
326 
327  LayerMultiplicityLayer0->Draw();
328  LayerMultiplicityLayer1->Draw("same");
329  LayerMultiplicityLayer2->Draw("same");
330 
331  LayerMultiplicityLayer0->SetTitle(";Chip multiplicity [Pixel];Count");
332 
333  TLegend *leg = new TLegend(.45, .55, .93, .93);
334  leg->AddEntry("", "#it{#bf{sPHENIX}} Simulation", "");
335  leg->AddEntry("", description, "");
336  leg->AddEntry("", Form("#mu_{MB} = %.1f, N_{Trig} = %d, #mu_{Noise} = %.1f", mu_MB, n_Trigger, mu_Noise), "");
337  leg->AddEntry(LayerMultiplicityLayer0, Form("MVTX Layer0 Chip4, <hit> = %.1f", LayerMultiplicityLayer0->GetMean()), "l");
338  leg->AddEntry(LayerMultiplicityLayer1, Form("MVTX Layer1 Chip4, <hit> = %.1f", LayerMultiplicityLayer1->GetMean()), "l");
339  leg->AddEntry(LayerMultiplicityLayer2, Form("MVTX Layer2 Chip4, <hit> = %.1f", LayerMultiplicityLayer2->GetMean()), "l");
340  leg->Draw();
341 
342  p = (TPad *) c1->cd(idx++);
343  c1->Update();
344  p->SetLogy();
345 
346  p->DrawFrame(0, 1. / NSample, MaxMultiplicity, 1)->SetTitle(";Chip multiplicity [Pixel];CCDF");
347  TLine *line = new TLine(0, 1e-3, MaxMultiplicity, 1e-3);
348  line->SetLineColor(kGray);
349  line->SetLineWidth(3);
350  line->Draw("same");
351 
352  MakeCDF(LayerMultiplicityLayer0)->Draw("same");
353  MakeCDF(LayerMultiplicityLayer1)->Draw("same");
354  MakeCDF(LayerMultiplicityLayer2)->Draw("same");
355 
356  SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), false);
357 
358  return cm;
359 }
360 
361 void OccupancySim( //
362  // const TString infile = "/sphenix/user/jinhuang/HF-jet/MVTX_Multiplicity/pp200MB_30cmVZ_Iter5/pp200MB_30cmVZ_Iter5_SUM.cfg_HFMLTriggerOccupancy.root", //
363  // const TString infile_trigger = "/sphenix/user/jinhuang/HF-jet/MVTX_Multiplicity/pp200MB_30cmVZ_Iter5/pp200MB_30cmVZ_Iter5_SUM.cfg_HFMLTriggerOccupancy.root", //
364  // const TString disc = "p+p MB, #sqrt{s} = 200 GeV" //
365  const TString infile = "/sphenix/user/jinhuang/HF-jet/MVTX_Multiplicity/AuAu200MB_30cmVZ_Iter5/AuAu200MB_30cmVZ_Iter5_SUM.xml_HFMLTriggerOccupancy.root", //
366  const TString infile_trigger = "/sphenix/user/jinhuang/HF-jet/MVTX_Multiplicity/AuAu200MB_10cmVZ_Iter5/AuAu200MB_10cmVZ_Iter5_SUM.xml_HFMLTriggerOccupancy.root", //
367  const TString disc = "Au+Au MB, #sqrt{s_{NN}} = 200 GeV" //
368 )
369 {
370  SetsPhenixStyle();
371  TVirtualFitter::SetDefaultFitter("Minuit2");
372 
373  description = disc;
374  _file0 = new TFile(infile);
375  assert(_file0->IsOpen());
376  _file_trigger = new TFile(infile_trigger);
377  assert(_file_trigger->IsOpen());
378 
379  _file0->cd();
380  Check(_file0);
381  chipMultiplicitySet_vec chipMultiplicity = MakeChipMultiplicitySet(_file0);
382 
383  _file_trigger->cd();
384  Check(_file_trigger);
385  chipMultiplicitySet_vec chipMultiplicityTrigger = MakeChipMultiplicitySet(_file_trigger);
386 
387  // AuAu
388  TriggerMultiplicity(chipMultiplicity, chipMultiplicityTrigger,
389  200e3 * 10e-6, 1, 1024 * 512 * 1e-6 + 2);
390 
391  TriggerMultiplicity(chipMultiplicity, chipMultiplicityTrigger,
392  200e3 * 15e-6, 0, 1024 * 512 * 1e-6 + 2);
393 
394  TriggerMultiplicity(chipMultiplicity, chipMultiplicityTrigger,
395  200e3 * 10e-6, 0, 1024 * 512 * 1e-6 + 2);
396 
397  TriggerMultiplicity(chipMultiplicity, chipMultiplicityTrigger,
398  200e3 * 5e-6, 0, 1024 * 512 * 1e-6 + 2);
399 
400  TriggerMultiplicity(chipMultiplicity, chipMultiplicityTrigger,
401  100e3 * 5e-6, 0, 1024 * 512 * 1e-6 + 2);
402 
403  TriggerMultiplicity(chipMultiplicity, chipMultiplicityTrigger,
404  50e3 * 5e-6, 0, 1024 * 512 * 1e-6 + 2);
405 
406  TriggerMultiplicity(chipMultiplicity, chipMultiplicityTrigger,
407  200e3 * 15e-6, 0, 1024 * 512 * 1e-5 + 2);
408 
409  TriggerMultiplicity(chipMultiplicity, chipMultiplicityTrigger,
410  200e3 * 10e-6, 0, 1024 * 512 * 1e-5 + 2);
411 
412  TriggerMultiplicity(chipMultiplicity, chipMultiplicityTrigger,
413  200e3 * 5e-6, 0, 1024 * 512 * 1e-5 + 2);
414 
415  // // pp
416  // TriggerMultiplicity(chipMultiplicity, chipMultiplicityTrigger,
417  // 13e6 * 10e-6, 1, 1024 * 512 * 1e-6);
418  //
419  // TriggerMultiplicity(chipMultiplicity, chipMultiplicityTrigger,
420  // 13e6 * 15e-6, 0, 1024 * 512 * 1e-6);
421  // TriggerMultiplicity(chipMultiplicity, chipMultiplicityTrigger,
422  // 13e6 * 10e-6, 0, 1024 * 512 * 1e-6);
423  // TriggerMultiplicity(chipMultiplicity, chipMultiplicityTrigger,
424  // 13e6 * 5e-6, 0, 1024 * 512 * 1e-6);
425  //
426  // TriggerMultiplicity(chipMultiplicity, chipMultiplicityTrigger,
427  // 13e6 * 15e-6, 0, 1024 * 512 * 1e-5);
428  // TriggerMultiplicity(chipMultiplicity, chipMultiplicityTrigger,
429  // 13e6 * 10e-6, 0, 1024 * 512 * 1e-5);
430  // TriggerMultiplicity(chipMultiplicity, chipMultiplicityTrigger,
431  // 13e6 * 5e-6, 0, 1024 * 512 * 1e-5);
432  // TriggerMultiplicity(chipMultiplicity, chipMultiplicityTrigger,
433  // 4e6 * 5e-6, 0, 1024 * 512 * 1e-5);
434 
435  //
436 }