Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
eicsmear_check_calorimeter.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file eicsmear_check_calorimeter.C
1 int
2 eicsmear_check_calorimeter( TString filename_mc,
3  TString filename_mc_smeared = "")
4 {
5  /* Uncomment this line when running without compilation */
6  gSystem->Load("libeicsmear");
7 
8  /* Open input files. */
9  TFile *file_mc = new TFile(filename_mc, "OPEN");
10  TFile *file_mc_smeared = new TFile(filename_mc_smeared, "OPEN");
11 
12  /* Get trees from files. */
13  TTree *tree = (TTree*)file_mc->Get("EICTree");
14  TTree *tree_smeared = (TTree*)file_mc_smeared->Get("Smeared");
15 
16  /* Add friend to match branches in trees. */
17  tree->AddFriend(tree_smeared);
18 
19  erhic::EventPythia *event = NULL;
20  Smear::Event *eventS = NULL;
21 
22  tree->SetBranchAddress("event", &event);
23  tree->SetBranchAddress("eventS", &eventS);
24 
25  /* Define base cut for reconstructed final state particles */
26  TCut cut_base("particles.KS==1 && Smeared.particles.p > 0 && (particles.id==11 || particles.id==22)");
27 
28  /* colors array */
29  unsigned colors[7] = {1,2,3,4,6,7,14};
30 
31  /* select how many eta settings to plot */
32  unsigned etas_plotmax = 7;
33 
34  /* Create vector of theta values to include */
35  vector< double > etas;
36  etas.push_back(-2.75);
37  etas.push_back(-2.25);
38  etas.push_back(-1.75);
39  etas.push_back(-0.25);
40  etas.push_back( 0.25);
41  etas.push_back( 1.75);
42  etas.push_back( 2.25);
43 
44  /* Create fit function */
45  TF1* f_momres = new TF1("f_momres", "sqrt( ([0])**2 + ([1]/(x**0.5))**2 )" );
46 
47  cout << "\nFit function: " << f_momres->GetTitle() << "\n" << endl;
48 
49  /* Create scratch canvas */
50  TCanvas *cscratch = new TCanvas("cscratch");
51 
52  /* Create framehistogram */
53  TH1F* hframe = new TH1F("hframe","",100,0,40);
54  hframe->GetYaxis()->SetRangeUser(0,0.15);
55  hframe->GetXaxis()->SetTitle("Energy (GeV/c)");
56  hframe->GetYaxis()->SetTitle("#sigma_{E}/E");
57 
58  /* Create 2Dhistogram */
59  unsigned nbinsp = 40;
60  TH2F* h2tmp = new TH2F("h2tmp","",nbinsp,0,(float)nbinsp,50,-1,1);
61  h2tmp->GetXaxis()->SetTitle("Energy (GeV)");
62  h2tmp->GetYaxis()->SetTitle("(#Delta E)/E");
63 
64  /* create combined canvas plot */
65  TCanvas *c1 = new TCanvas();
66  hframe->Draw();
67 
68  /* Create legend */
69  TLegend* leg_eta = new TLegend( 0.25, 0.40, 0.45, 0.90);
70 
71  /* Create resolution-vs-momentum plot with fits for each selected theta value */
72  for ( int i = 0; i < etas.size(); i++ )
73  {
74  /* Switch to scratch canvas */
75  cscratch->cd();
76 
77  double eta = etas.at(i);
78 
79  cout << "\n***Eta = " << eta << endl;
80 
81  /* Define range of theta because float comparison with fixed value doesn't work
82  too well for cuts in ROOT trees */
83  double eta_min = eta - 0.25;
84  double eta_max = eta + 0.25;
85 
86  /* Cut for tree */
87  TCut cutx( Form("particles.p > 1 && ( particles.eta > %f && particles.eta < %f )", eta_min, eta_max) );
88  cout << cutx.GetTitle() << endl;
89 
90  /* Fill 2D histogram and fit slices with gaussian to get track resolution */
91  h2tmp->Reset();
92  tree->Draw("(Smeared.particles.E-particles.E)/particles.E:particles.E >> h2tmp",cut_base && cutx,"");
93  h2tmp->FitSlicesY();
94 
95  /* Get histogram with sigma from gauss fits */
96  TH1D* h2tmp_2 = (TH1D*)gDirectory->Get("h2tmp_2");
97 
98  /* Fill fit parameters in TGraphErrors */
99  TGraphErrors *gres = new TGraphErrors(nbinsp);
100  gres->SetMarkerColor(colors[i]);
101  for ( unsigned bini = 1; bini < nbinsp; bini++ )
102  {
103  double sigm_i = h2tmp_2->GetBinContent(bini);
104  double sigm_err_i = h2tmp_2->GetBinError(bini);
105 
106  gres->SetPoint( bini-1, h2tmp->GetXaxis()->GetBinCenter(bini), sigm_i );
107  gres->SetPointError( bini-1, 0, sigm_err_i );
108  }
109 
110  /* Only plot first few lines; if not plotting, still do the fit */
111  if ( i < etas_plotmax )
112  {
113  /* Add graph to legend */
114  leg_eta->AddEntry(gres, Form("#eta = %.1f", eta), "P");
115 
116  /* Add graph to plot */
117  c1->cd();
118  gres->Draw("Psame");
119  f_momres->SetLineColor(colors[i]);
120  gres->Fit(f_momres);
121  }
122  else
123  {
124  gres->Fit(f_momres);
125  }
126  }
127 
128  /* Draw legend */
129  TCanvas *c2 = new TCanvas();
130  leg_eta->Draw();
131 
132  return 0;
133 }