StRoot  1
GetCentrality.C
1 #include <iostream>
2 #include "TMath.h"
3 #include "TRandom3.h"
4 #include "TCanvas.h"
5 #include "TF1.h"
6 #include "TH1D.h"
7 #include "TH1D.h"
8 #include "TGraph.h"
9 #include "TAxis.h"
10 #include "TGraph2D.h"
11 #include "TGraphErrors.h"
12 #include "TH1F.h"
13 #include "TFile.h"
14 #include "TNtuple.h"
15 #include <numeric>
16 #include <sstream>
17
18 void GetCentrality() {
19  TFile *f0 = new TFile("glauberfit.root");
20  TH1D *sim = (TH1D *)f0->Get("hRefMultSim");
21  TH1D *data = (TH1D *)f0->Get("hRefMultCorr");
22  //sim->Scale(1.0/sim->Integral(50,500));
23  //data->Scale(1.0/data->Integral(50,500));
24
25
26  int centralitybin[16][2];
27  double integral = sim->Integral();
28  for(int cent=0; cent<16; cent++){
29  double distance = 1000.0; //Initialize the distance from 5% cut to be large
30  double fraction = 0.05*((double)cent+1.0);
31  //For the most central bins, integrate the data
32  if(cent<4){
33  for(int i=0; i<data->GetNbinsX(); i++){
34  double thisfraction = (data->Integral(i,500))/(integral);
35  double thisdistance = TMath::Abs(thisfraction - fraction);//This is how far the current integral
36  // is from the desired fraction.
37  if(thisdistance>distance){ //If the distance from the fraction increased, then the previous bin was your desired cut
38  if(cent==0){centralitybin[0][0]=500; centralitybin[0][1]=i-1;centralitybin[1][0]=i-2;}
39  else{centralitybin[cent][1]=i-1;centralitybin[cent+1][0]=i-2;}
40  break;
41  }
42  else{distance=thisdistance;}
43  }
44  }
45  //For more peripheral bins, integrate the simulated distribution
46  else{
47  int newmaxbin = centralitybin[3][1]-1;//Integrate the glauber up to this bin
48  double zeroToTwentyIntegral = data->Integral(newmaxbin+1,500);
49  for(int i=0; i<newmaxbin; i++){
50  double thisfraction=(sim->Integral(i,newmaxbin)+zeroToTwentyIntegral)/(integral);
51  double thisdistance = TMath::Abs(thisfraction - fraction);
52  if(thisdistance>distance){
53  if(cent==15){centralitybin[15][1]=i-1;}
54  else{centralitybin[cent][1]=i-1;centralitybin[cent+1][0]=i-2;}
55  break;
56  }
57  else{distance=thisdistance;}
58  }
59  }
60  }
61
62  //Print out all useful imformation
63  cout<<"******* 16 Bins *******"<<endl;
64  cout<<"High bins"<<endl;
65  for(int i=0; i<16; i++){
66  cout<<data->GetXaxis()->GetBinCenter(centralitybin[i][0])+0.5<<endl;
67  }
68  cout<<"Low bins"<<endl;
69  for(int i=0; i<16; i++){
70  cout<<data->GetXaxis()->GetBinCenter(centralitybin[i][1])-0.5<<endl;
71  }
72  cout<<"Integrals"<<endl;
73  for(int i=0; i<16; i++){
74  if(i<4)cout<<(data->Integral(centralitybin[i][1],centralitybin[i][0]))/(sim->Integral())<<endl;
75  else cout<<(sim->Integral(centralitybin[i][1],centralitybin[i][0]))/(sim->Integral())<<endl;
76  }
77  cout<<"Cumulative Integrals"<<endl;
78  for(int i=0; i<16; i++){
79  if(i<4)cout<<(data->Integral(centralitybin[i][1],500))/(sim->Integral())<<endl;
80  else cout<<(sim->Integral(centralitybin[i][1],centralitybin[3][1]-1)+data->Integral(centralitybin[3][1],500))/(sim->Integral())<<endl;
81  }
82  cout<<"Efficiencies"<<endl;
83  for(int i=0; i<16; i++){
84  double thiseff = (data->Integral(centralitybin[i][1],centralitybin[i][0]))/(sim->Integral(centralitybin[i][1],centralitybin[i][0]));
85  if(thiseff>1.0) thiseff=1.0;
86  if(i<4)cout<<1<<endl;
87  else cout<<thiseff<<endl;
88  }
89  /*
90  cout<<"<refMultCorr>"<<endl;
91  for(int i=0; i<16; i++){
92  double thisavg = data->GetBinContent(centralitybin[i][1],centralitybin[i][0]))/(sim->Integral(centralitybin[i][1],centralitybin[i][0]));
93  if(thiseff>1.0) thiseff=1.0;
94  if(i<4)cout<<1<<endl;
95  else cout<<thiseff<<endl;
96  }
97  */
98  cout<<"******* 9 Bins *******"<<endl;
99  cout<<"High bins"<<endl;
100  for(int i=0; i<9; i++){
101  if(i==0 || i==1){cout<<data->GetXaxis()->GetBinCenter(centralitybin[i][0])+0.5<<endl;}
102  else{cout<<data->GetXaxis()->GetBinCenter(centralitybin[2*i-1][0])+0.5<<endl;}
103  }
104  cout<<"Low bins"<<endl;
105  for(int i=0; i<9; i++){
106  if(i==0 || i==1){cout<<data->GetXaxis()->GetBinCenter(centralitybin[i][1])-0.5<<endl;}
107  else{cout<<data->GetXaxis()->GetBinCenter(centralitybin[2*i-1][1])-0.5<<endl;}
108  }
109 }