StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
myProjectProfile.C
1 #include "TH1.h"
2 #include "TProfile.h"
3 #include "TProfile2D.h"
4 #include "iostream.h"
5 
6 //absoluteBould =1 make, one can use it like 2.6<|eta|<4
7 //reversedSign =1 make the addition of profile as ( negProf *-1. + posProf)
8 
9 
10 TH1D* myProjectionX(TProfile2D* prof2d, double lowBound, double upBound, bool absoluteBound=0, bool reversedSign=0 ){
11 
12 
13 
14 
15  TString st(prof2d->GetName());
16  st +="_X";
17 
18 
19  TH1D* projx = prof2d->ProjectionX(st.Data(), prof2d->GetYaxis()->FindBin(lowBound),prof2d->GetYaxis()->FindBin(upBound));
20  projx->Reset();
21 
22  TH1D* projy = prof2d->ProjectionY();
23 
24 
25 
26 
27  for (int i = 0; i<projx->GetNbinsX(); i++) {
28 
29  double contentSum=0.;
30  double entriesSum=0.;
31  double contentSqrSum=0.;
32 
33 
34  for (int j=0; j<projy->GetNbinsX(); j++) {
35 
36 
37  float BinHighEdge = (projy->GetBinLowEdge(j) + projy->GetBinWidth(j));
38  float BinLowEdge = projy->GetBinLowEdge(j);
39 
40  if (absoluteBound==0){
41  if ( BinHighEdge < lowBound) continue;
42  if ( BinLowEdge > upBound) continue;
43  } else if (absoluteBound){
44 
45  if ( BinHighEdge < 0. && BinLowEdge < 0. && BinHighEdge < -upBound) continue;
46  if ( BinHighEdge < 0. && BinLowEdge < 0. && BinLowEdge > -lowBound) continue;
47  if ( BinHighEdge > 0. && BinLowEdge > 0. && BinHighEdge < lowBound) continue;
48  if ( BinHighEdge > 0. && BinLowEdge > 0. && BinLowEdge > upBound) continue;
49  if ( BinHighEdge > 0. && BinLowEdge < 0. && (BinHighEdge-BinLowEdge)<lowBound) continue;
50  }
51 
52 
53 
54 
55 
56  double entries = prof2d->GetBinEntries(prof2d->GetBin(i,j));
57  if (entries==0) continue;
58  double content = prof2d->GetBinContent(i,j);
59  if (reversedSign && (BinHighEdge+BinLowEdge)/2. < 0.) content *=-1.;
60  double error = prof2d->GetBinError(i,j);
61 
62  entriesSum += entries;
63  contentSum += content * entries;
64 
65  if (entries > 1) {
66  contentSqrSum += entries * (entries * error * error + content * content);
67  } else if (entries ==1) { //special case to calculate the correct error
68 
69  contentSqrSum += content*content;
70  }
71  }
72 
73  if (entriesSum) {
74 
75  projx->SetBinContent(i,contentSum/entriesSum);
76  projx->SetBinError(i,(sqrt((contentSqrSum/entriesSum) - pow ( ( contentSum/entriesSum),2.)) / sqrt(entriesSum)));
77 
78  }
79 
80  }
81 
82  if (projy) delete projy;
83  return projx;
84 
85 }
86 
87 TH1D* myProjectionY(TProfile2D* prof2d, double lowBound, double upBound ,bool absoluteBound=0, bool reversedSign=0 ){
88 
89  TString st(prof2d->GetName());
90  st +="_Y";
91 
92 
93 
94  TH1D* projy = prof2d->ProjectionY(st.Data(), prof2d->GetXaxis()->FindBin(lowBound),prof2d->GetXaxis()->FindBin(upBound));
95  projy->Reset();
96 
97  TH1D* projx = prof2d->ProjectionX();
98 
99 
100  for (int j=0; j<projy->GetNbinsX(); j++) {
101 
102  double contentSum=0.;
103  double entriesSum=0.;
104  double contentSqrSum=0.;
105 
106  for (int i = 0; i<projx->GetNbinsX(); i++) {
107 
108 
109  float BinHighEdge = (projx->GetBinLowEdge(i) + projx->GetBinWidth(i));
110  float BinLowEdge = projx->GetBinLowEdge(i);
111 
112  if (absoluteBound==0){
113  if ( BinHighEdge < lowBound) continue;
114  if ( BinLowEdge > upBound) continue;
115  } else if (absoluteBound){
116 
117  if ( BinHighEdge < 0. && BinLowEdge < 0. && BinHighEdge < -upBound) continue;
118  if ( BinHighEdge < 0. && BinLowEdge < 0. && BinLowEdge > -lowBound) continue;
119  if ( BinHighEdge > 0. && BinLowEdge > 0. && BinHighEdge < lowBound) continue;
120  if ( BinHighEdge > 0. && BinLowEdge > 0. && BinLowEdge > upBound) continue;
121  if ( BinHighEdge > 0. && BinLowEdge < 0. && (BinHighEdge-BinLowEdge)<lowBound) continue;
122  }
123 
124 
125 
126 
127  double entries = prof2d->GetBinEntries(prof2d->GetBin(i,j));
128  if (entries==0) continue;
129  double content = prof2d->GetBinContent(i,j);
130  if (reversedSign && (BinHighEdge+BinLowEdge)/2. < 0.) content *=-1.; // for directed flow
131  double error = prof2d->GetBinError(i,j);
132 
133  entriesSum += entries;
134  contentSum += content * entries;
135 
136  if (entries > 1) {
137  contentSqrSum += entries * (entries * error * error + content * content);
138  } else if (entries ==1) { //special case to calculate the correct error
139 
140  contentSqrSum += content*content;
141  }
142  }
143 
144  if (entriesSum) {
145 
146  projy->SetBinContent(j,contentSum/entriesSum);
147  projy->SetBinError(j,(sqrt((contentSqrSum/entriesSum) - pow ( ( contentSum/entriesSum),2.)) / sqrt(entriesSum)));
148 
149  }
150 
151  }
152  if (projx) delete projx;
153  return projy;
154 
155 }