StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
minBias.C
1 //
3 // $Id: minBias.C,v 1.15 2007/02/06 19:00:51 posk Exp $
4 //
5 // Author: Art Poskanzer and Alexander Wetzler, Mar 2001
6 // Kirill Filimonov treated the one count case
7 // Description: Macro to add histograms together.
8 // The v histograms will be added with yield weighting.
9 // First file anaXX.root given by first run number XX.
10 // Last file anaXX.root given by last run number XX.
11 // Output file given by output run number.
12 // For 10-50% centrality: minBias(X4,X7,X0)
13 //
14 //
16 #ifndef __CINT__
17 #include <fstream.h>
18 #include "TSystem.h"
19 #include <TFile.h>
20 #include "TH1.h"
21 #include "TH2.h"
22 #include "TProfile.h"
23 #include "TKey.h"
24 #include "TObject.h"
25 #endif
26 
27 void minBias(Int_t firstRunNo, Int_t lastRunNo, Int_t outputRunNo=99) {
28 
29  const int nCens = lastRunNo - firstRunNo + 1;
30  int nSels = 2;
31  const int nHars = 4;
32  char fileName[80];
33  TFile* histFile[nCens+1];
34  TH1* hist[nCens+1];
35  TH2* yieldPartHist[nCens];
36  TH1* yieldPartHistPt[nCens];
37  TH1* yieldPartHistEta[nCens];
38  Float_t yieldTotal[nCens];
39 
40  // names of histograms to be added with weighting
41  const char* baseName[] = {
42  "Flow_Res_",
43  "Flow_v2D_", // must be before vEta and vPt
44  "Flow_vEta_",
45  "Flow_vPt_",
46  "Flow_v_",
47  "FlowLYZ_r0theta_",
48  "FlowLYZ_Vtheta_",
49  "FlowLYZ_V_",
50  "FlowLYZ_vr0_",
51  "FlowLYZ_vEta_",
52  "FlowLYZ_vPt_",
53  "FlowLYZ_v_",
54  "Flow_v2D_ScalarProd_",
55  "Flow_vEta_ScalarProd_",
56  "Flow_vPt_ScalarProd_",
57  "Flow_v_ScalarProd_",
58  "Flow_Cumul_vEta_Order2_",
59  "Flow_Cumul_vPt_Order2_",
60  "Flow_Cumul_v_Order2_",
61  "Flow_Cumul_vEta_Order4_",
62  "Flow_Cumul_vPt_Order4_",
63  "Flow_Cumul_v_Order4_"
64  };
65  const int nNames = sizeof(baseName) / sizeof(char*);
66 
67  // open the files
68  for (int n = 0; n < nCens; n++) {
69  sprintf(fileName, "ana%2d.root", firstRunNo + n);
70  cout << " file name = " << fileName << endl;
71  histFile[n] = new TFile(fileName);
72  if (!histFile[n]) {
73  cout << "### Can't find file " << fileName << endl;
74  return;
75  }
76  }
77  sprintf(fileName, "ana%2d.root", outputRunNo);
78  cout << " output file name = " << fileName << endl << endl;
79  histFile[nCens] = new TFile(fileName, "RECREATE");
80 
81  // add all histograms with no weighting
82  TKey* key;
83  TObject* obj;
84  TIter nextkey(histFile[0]->GetListOfKeys());
85  while (key = (TKey*)nextkey()) {
86  histFile[0]->cd();
87  //key->ls();
88  obj = key->ReadObj();
89  const char* objName;
90  if (obj->InheritsFrom("TH1")) { // TH1 or TProfile
91  hist[0] = (TH1*)obj;
92  objName = key->GetName();
93  cout << "hist name= " << objName << endl;
94  for (int n = 1; n < nCens; n++) {
95  hist[1] = (TH1*)histFile[n]->Get(objName);
96  hist[0]->Add(hist[1]);
97  }
98  }
99  histFile[nCens]->cd();
100  obj->Write(objName);
101  delete obj;
102  obj = NULL;
103  }
104 
105  // get yield histogram
106  cout<<endl;
107  for (int n = 0; n < nCens; n++) {
108  yieldPartHist[n] = dynamic_cast<TH2*>(histFile[n]->Get("Flow_YieldPart2D"));
109  yieldPartHistPt[n] = dynamic_cast<TH1*>(histFile[n]->Get("FlowLYZ_YieldPartPt"));
110  yieldPartHistEta[n] = dynamic_cast<TH1*>(histFile[n]->Get("FlowLYZ_YieldPartEta"));
111  if (yieldPartHistPt[n]) { yieldTotal[n] = yieldPartHistPt[n]->Integral(); }
112  }
113 
114  for (int pageNumber = 0; pageNumber < nNames; pageNumber++ ) {
115  bool twoD = kFALSE;
116  if (strstr(baseName[pageNumber],"v2D")) twoD = kTRUE;
117 
118  for (int selN = 0; selN < nSels; selN++) {
119 
120  // no harmonics
121  bool noHars = kFALSE;
122  int nHar = nHars;
123  if (strstr(baseName[pageNumber],"_v_")!=0 ||
124  strstr(baseName[pageNumber],"Res")!=0 ||
125  strstr(baseName[pageNumber],"_V_")!=0 ||
126  strstr(baseName[pageNumber],"_vr0_")!=0) {
127  noHars = kTRUE;
128  nHar = 1;
129  }
130  for (int harN = 0; harN < nHar; harN++) {
131 
132  // construct histName
133  TString* histName = new TString(baseName[pageNumber]);
134  *histName += "Sel";
135  *histName += selN+1;
136  if (!noHars) {
137  *histName += "_Har";
138  *histName += harN+1;
139  }
140 
141  // get the histograms
142  for (int n = 0; n < nCens+1; n++) {
143  hist[n] = dynamic_cast<TH1*>(histFile[n]->Get(histName->Data()));
144  }
145  if (hist[0]) { cout << "hist name= " << histName->Data() << endl; }
146 
147  const int lastHist = nCens;
148  int nBins; // set by 2D of centrality lastHist
149  if (!hist[lastHist]) continue;
150  int xBins = hist[lastHist]->GetNbinsX();
151  int yBins;
152  TH1F *yieldY, *yieldPt;
153  if (twoD) {
154  yBins = hist[lastHist]->GetNbinsY();
155  nBins = xBins + (xBins + 2) * yBins;
156  float yMax = hist[lastHist]->GetXaxis()->GetXmax();
157  float yMin = hist[lastHist]->GetXaxis()->GetXmin();
158  yieldY = new TH1F("Yield_Y", "Yield_Y", xBins, yMin, yMax);
159  yieldY->SetXTitle("Rapidity");
160  yieldY->SetYTitle("Counts");
161  float ptMax = hist[lastHist]->GetYaxis()->GetXmax();
162  yieldPt = new TH1F("Yield_Pt", "Yield_Pt", yBins, 0., ptMax);
163  yieldPt->SetXTitle("Pt (GeV/c)");
164  yieldPt->SetYTitle("Counts");
165  } else {
166  nBins = xBins + 2;
167  }
168 
169  // loop over the bins
170  if (strstr(baseName[pageNumber],"Res") ||
171  strstr(baseName[pageNumber],"Cumul")) { // with error weighting
172  cout << " With error weighting" << endl;
173  float content, error, errorSq, meanContent, meanError, weight;
174  for (int bin = 0; bin < nBins; bin++) {
175  meanContent = 0.;
176  meanError = 0.;
177  weight = 0.;
178  for (int n = 0; n < nCens; n++) {
179  if (hist[n]) {
180  content = hist[n]->GetBinContent(bin);
181  error = hist[n]->GetBinError(bin);
182  errorSq = error * error;
183  if (errorSq > 0.) {
184  meanContent += content / errorSq;
185  weight += 1. / errorSq;
186  }
187  }
188  }
189  if (weight > 0.) {
190  meanContent /= weight;
191  meanError = sqrt(1. / weight);
192  hist[nCens]->SetBinContent(bin, meanContent);
193  hist[nCens]->SetBinError(bin, meanError);
194  }
195  }
196  } else { // with yield weighting
197  cout << " With yield weighting" << endl;
198  float v, verr, content, error, yield, y, pt;
199  double vSum, vSum2, error2sum, yieldSum, vRms, verrRms, vSumRms, yieldSumRms, vSumRms2;
200  for (int bin = 0; bin < nBins; bin++) {
201  v = 0.;
202  verr = 0.;
203  vSum = 0.;
204  vSum2 = 0.;
205  content = 0.;
206  error = 0.;
207  error2sum = 0.;
208  yield = 0.;
209  yieldSum = 0.;
210  vRms = 0.;
211  verrRms = 0.;
212  vSumRms = 0.;
213  yieldSumRms = 0.;
214  vSumRms2 = 0.;
215  for (int n = 0; n < nCens; n++) {
216  if (hist[n]) {
217  if (strstr(baseName[pageNumber],"LYZ")) {
218  if (strstr(histName->Data(), "vEta")) {
219  yield = yieldPartHistEta[n]->GetBinContent(bin);
220  } else if (strstr(histName->Data(), "vPt")) {
221  yield = yieldPartHistPt[n]->GetBinContent(bin);
222  } else { // r0theta, Vtheta, _V_, vr0, _v_, Gtheta
223  yield = yieldTotal[n];
224  }
225  } else {
226  if (strstr(histName->Data(),"v2D")) {
227  yield = yieldPartHist[n]->GetBinContent(bin);
228  } else if (strstr(histName->Data(),"vEta")) {
229  yield = yieldPartHist[n]->Integral(bin, bin, 1, yBins);
230  if (selN==0 && harN==0) {
231  y = yieldPartHist[n]->GetXaxis()->GetBinCenter(bin);
232  yieldY->Fill(y, yield);
233  }
234  } else if (strstr(histName->Data(),"vPt")) {
235  yield = yieldPartHist[n]->Integral(1, xBins, bin, bin);
236  if (selN==0 && harN==0) {
237  pt = yieldPartHist[n]->GetYaxis()->GetBinCenter(bin);
238  yieldPt->Fill(pt, yield);
239  }
240  } else { // _v
241  yield = yieldPartHist[n]->Integral();
242  }
243  }
244  v = hist[n]->GetBinContent(bin);
245  if (yield==1) { // special case to calculate the correct error
246  vSumRms += v;
247  yieldSumRms += yield;
248  vSumRms2 += v*v;
249  } else {
250  verr = hist[n]->GetBinError(bin);
251  if (v != 0 ) {
252  yieldSum += yield;
253  vSum += yield * v;
254  vSum2 += v * v * yield;
255  error2sum += pow(yield * verr, 2.);
256  }
257  }
258  }
259  }
260 
261  if (yieldSumRms) vRms = vSumRms / yieldSumRms;
262  if (yieldSumRms>1) {
263  verrRms = sqrt(vSumRms2 - vSumRms*vSumRms / yieldSumRms)
264  / yieldSumRms;
265  yieldSum += yieldSumRms;
266  vSum += vSumRms;
267  error2sum += pow(yieldSumRms * verrRms, 2.);
268  vSum2 += vRms * vRms * yieldSumRms;
269  }
270 
271  if (yieldSum) {
272  content = vSum / yieldSum;
273  if (yieldSumRms==1) {
274  error = sqrt(error2sum + vSum2 - vSum*vSum / yieldSum) / yieldSum;
275  error = yieldSum / (yieldSum+1) * sqrt(error*error +
276  (vRms - content)*(vRms - content) / (yieldSum * (yieldSum+1)));
277  } else {
278  error = sqrt(error2sum + vSum2 - vSum*vSum/yieldSum) / yieldSum;
279  }
280  hist[nCens]->SetBinContent(bin, content);
281  hist[nCens]->SetBinError(bin, error);
282  }
283  }
284  }
285  delete histName;
286  }
287  }
288  }
289 
290  //histFile[nCens]->ls();
291  histFile[nCens]->Write(0, TObject::kOverwrite);
292  histFile[nCens]->Close();
293  delete histFile[nCens];
294 
295 }
296 
298 //
299 // $Log: minBias.C,v $
300 // Revision 1.15 2007/02/06 19:00:51 posk
301 // In Lee Yang Zeros method, introduced recentering of Q vector.
302 // Reactivated eta symmetry cut.
303 //
304 // Revision 1.13 2006/03/22 22:02:07 posk
305 // Updates to macros.
306 //
307 // Revision 1.12 2006/02/24 18:36:04 posk
308 // Updated for LeeYangZerosMaker histograms.
309 //
310 // Revision 1.11 2004/03/01 22:43:42 posk
311 // Changed some "->" to ".".
312 //
313 // Revision 1.10 2003/08/26 21:10:12 posk
314 // Calculates v8 if nHars=8.
315 //
316 // Revision 1.9 2003/06/27 21:25:44 posk
317 // v4 and v6 are with repect to the 2nd harmonic event plane.
318 //
319 // Revision 1.8 2003/03/11 23:03:07 posk
320 // Includes scalar product and cumulant hists.
321 //
322 // Revision 1.7 2002/06/11 21:54:15 posk
323 // Kirill's further correction to minBias.C for bins with one count.
324 //
325 // Revision 1.6 2002/05/21 18:42:18 posk
326 // Kirill's correction to minBias.C for bins with one count.
327 //
328 // Revision 1.5 2001/11/09 21:15:00 posk
329 // Switched from CERNLIB to TMath. Using global dca instead of dca.
330 //
331 // Revision 1.4 2001/05/22 20:05:47 posk
332 // Now outputs a hist.root file.
333 // The v values are averaged with yield weighting.
334 //
335 // Revision 1.3 2000/09/29 22:53:17 posk
336 // More histograms.
337 //
338 // Revision 1.2 2000/09/26 20:54:11 posk
339 // Updated documentation.
340 //
341 // Revision 1.1 2000/09/26 00:19:43 posk
342 // New macro to add centrality-selected histograms with proper weights, to make
343 // minimum bias histogram.
344 //
345 //
Generated by   doxygen 1.8.5