StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
doFlowSumFirstPass.C
1 //
3 // $Id: doFlowSumFirstPass.C,v 1.2 2006/07/06 16:58:37 posk Exp $
4 //
5 // Makes temporary root.files.<cenNo> containing lists of flow.firstPassLYZNew.root files
6 // in all subdirectories of outDir.
7 // outDir is a link in this directory.
8 // The cenNo is the firstCenNo + cen.
9 // The link argument is the prefix of the subdirectory name.
10 // Adds the histograms together by yield weighted averaging for all files for this centrality.
11 // Writes out files called flow.firstPassLYZ<cenNo>.root.
12 // Copies them back to the subdirectories as flow.firstPassLYZNew.root.
13 // This macro must be in your directory containing the outDir link
14 // and run between the first and second passes.
15 //
16 // by Art Poskanzer, Feb 2006
17 //
19 
20 #include <iostream.h>
21 #include <iomanip.h>
22 #include <fstream.h>
23 #include "TH1.h"
24 #include "TFile.h"
25 #include "TComplex.h"
26 
27 void doFlowSumFirstPass(char* link, int firstCenNo) {
28 
29  int nCens = 9; // 9
30  Bool_t fromG = kTRUE; // the prefered setting
31 
32  char rootFileName[80];
33  char rootDirName[80];
34  char listFileName[80];
35  char rootOutName[80];
36  char lsCommand[80];
37  char firstPassCommand[80];
38  char rmCommand[80];
39  char cpCommand[80];
40  TFile* histFile[1000];
41  Float_t j01 = 2.405;
42  int maxTheta;
43  const int nSels = 2;
44  const int nHars = 2;
45 
46  // names of histograms to be averaged with weighting (no profiles and no 2D hists)
47  const char* baseName[] = {
48  "FlowLYZ_r0theta_",
49  "FlowLYZ_Vtheta_", // must come after r0theta
50  "FlowLYZ_Gtheta", // first G
51  "FlowReGtheta",
52  "FlowImGtheta" // after ReG
53  };
54  const int nNames = sizeof(baseName) / sizeof(char*);
55 
56  for (int cen = 0; cen < nCens; cen++) {
57  int cenNo = firstCenNo + cen;
58 
59  // open the input files
60  Int_t nFile = 0;
61  sprintf(listFileName, "root.files.%2d", cenNo);
62  sprintf(lsCommand, "ls -1 outDir/%s-*-%2d/flow.firstPassLYZNew.root > %s", link, cenNo,
63  listFileName);
64  //sprintf(lsCommand, "ls -1 flow.firstPassLYZNew.root > %s", listFileName);
65  system(lsCommand);
66  fstream ascii_in;
67  ascii_in.open(listFileName, ios::in);
68  ascii_in >> rootFileName;
69  while(!ascii_in.eof()) {
70  if (strstr(rootFileName,".root")==0) continue;
71  histFile[nFile] = new TFile(rootFileName);
72  //if (nFile==0) { histFile[nFile]->ls(); }
73  strcpy(rootDirName, rootFileName);
74  char* cp = strstr(rootDirName, "flow.");
75  if (cp) *cp = '\0'; // truncate to directory name
76  sprintf(firstPassCommand, "test -f %s", rootFileName);
77  if (system(firstPassCommand)) {
78  cout << "### No first pass for " << rootDirName << endl;
79  //delete histFile[nFile];
80  continue;
81  }
82  cout << nFile+1 << " directory name = " << rootDirName << endl;
83  ascii_in >> rootFileName;
84  nFile++;
85  }
86  const Int_t nFiles = nFile;
87  cout << "input files: " << nFiles << endl << endl;
88  if (nFiles == 0) { return; }
89 
90  // the following 3 lines will prevent compiling this macro
91  TH1D* hist[nFiles+1]; // last file is the output
92  TH1D* hist_r0theta[nSels][nHars];
93  TH1F* yieldMultHist[nFiles+1];
94  Float_t yieldTotal[nFiles];
95 
96  // open the output file
97  sprintf(rootOutName, "flow.firstPassLYZ%2d.root", cenNo);
98  cout << "tempory output file name = " << rootOutName << endl << endl ;
99  histFile[nFiles] = new TFile(rootOutName, "RECREATE"); // last file is new out file
100 
101  // get multiplicity histograms
102  for (int n = 0; n < nFiles; n++) {
103  yieldMultHist[n] = dynamic_cast<TH1F*>(histFile[n]->Get("FlowLYZ_Mult"));
104  if (!yieldMultHist[n]) {
105  cout << "### dynamic cast can't find histogram FlowLYZ_Mult" << endl;
106  return;
107  }
108  yieldTotal[n] = yieldMultHist[n]->GetEntries() * yieldMultHist[n]->GetMean();
109  cout << n+1 << ": " << setprecision(7) << yieldTotal[n] << " particles" << endl;
110  }
111 
112  // add mult hists together for output
113  histFile[nFiles]->cd();
114  yieldMultHist[nFiles] = (TH1F*)yieldMultHist[0]->Clone();
115  for (int n = 1; n < nFiles; n++) {
116  yieldMultHist[nFiles]->Add(yieldMultHist[n]);
117  }
118  cout << endl;
119 
120  // Get maxTheta
121  if (cen==0) {
122  TString histNameNtheta("FlowLYZ_r0theta_Sel1_Har2");
123  TH1D* histNtheta = (TH1D*)histFile[0]->Get(histNameNtheta);
124  if (!histNtheta) {
125  cout << "### Can't find file " << histNameNtheta << endl;
126  return;
127  }
128  maxTheta = histNtheta->GetNbinsX();
129  cout << "maxTheta = " << maxTheta << endl << endl;
130  const int thetas = maxTheta + 2;
131  }
132  TH1D* histReG[nSels][nHars][thetas];
133  TH1D* histG[nSels][nHars][thetas];
134 
135  //if (fromG) { cout << "r0:\t\t\t from V \t from G \t change" << endl; }
136  float r0V[nSels][nHars][thetas];
137  float Xlast, X0, Xnext;
138  double reG, imG, Glast, G0, Gnext, GnextNext;
139  TComplex Gtheta;
140 
141  // with yield weighted averaging
142  for (int pageNumber = 0; pageNumber < nNames; pageNumber++ ) {
143  for (int selN = 0; selN < nSels; selN++) {
144  for (int harN = 0; harN < nHars; harN++) {
145  for (int thetaN = 0; thetaN < maxTheta; thetaN++) {
146 
147  // construct histName
148  TString histName(baseName[pageNumber]);
149  if (strstr(histName.Data(), "G")) {
150  histName += thetaN;
151  histName += "_Sel";
152  histName += selN+1;
153  histName += "_Har";
154  histName += harN+1;
155  } else {
156  histName += "Sel";
157  histName += selN+1;
158  histName += "_Har";
159  histName += harN+1;
160  }
161  //cout << "hist name= " << histName << endl;
162 
163  // get the histograms
164  for (int n = 0; n < nFiles; n++) {
165  hist[n] = dynamic_cast<TH1D*>(histFile[n]->Get(histName));
166  if (!hist[n]) {
167  cout << "### dynamic cast can't find file " << n << ": histogram " <<
168  histName << endl;
169  return;
170  }
171  }
172 
173  // make the output hists
174  histFile[nFiles]->cd();
175  hist[nFiles] = new TH1D(*(hist[0]));
176 
177  // get pointer to the r0 output hist
178  if (strstr(histName.Data(), "r0theta")) {
179  hist_r0theta[selN][harN] = dynamic_cast<TH1D*>(histFile[nFiles]->Get(histName));
180  if (!hist_r0theta[selN][harN]) {
181  cout << "### dynamic castan't find hist " << histName << endl;
182  return;
183  }
184  }
185 
186  // get pointer to the G output hist
187  if (strstr(histName.Data(), "_G")) {
188  histG[selN][harN][thetaN] = dynamic_cast<TH1D*>(histFile[nFiles]->Get(histName));
189  if (!histG[selN][harN][thetaN]) {
190  cout << "### dynamic cast can't find hist " << histName << endl;
191  return;
192  }
193  }
194 
195  // get the number of bins for this hist
196  if (!hist[0]) continue;
197  int nBins = hist[0]->GetNbinsX() + 2;
198 
199  // loop over the bins
200  float yield, r0fromV;
201  double content, error, Vtheta, VthetaErr, r0, r0VErr;
202  double v, vSum, v2Sum, error2Sum, yieldSum, yieldSum2, yield2Sum;
203  for (int bin = 0; bin < nBins; bin++) {
204  v = 0.;
205  vSum = 0.;
206  v2Sum = 0.;
207  content = 0.;
208  error = 0.;
209  error2Sum = 0.;
210  yield = 0.;
211  yieldSum = 0.;
212  yield2Sum = 0.;
213  for (int n = 0; n < nFiles; n++) { // loop over the files
214  if (hist[n]) {
215  yield = yieldTotal[n];
216  v = hist[n]->GetBinContent(bin);
217  if (v != 0. && yield > 1.) { // error is wrong for one count
218  yieldSum += yield;
219  yield2Sum += yield*yield;
220  vSum += yield * v;
221  v2Sum += yield*yield * v*v;
222  error2Sum += pow(yield * hist[n]->GetBinError(bin), 2.);
223  }
224  }
225  } // nFiles
226  if (yieldSum) {
227  content = vSum / yieldSum;
228  yieldSum2 = yieldSum*yieldSum;
229  v2Sum /= yieldSum2;
230  yield2Sum /= yieldSum2;
231  error2Sum /= yieldSum2;
232  if (strstr(histName.Data(), "G")) {
233  error = 0.;
234  } else {
235  error = sqrt(error2Sum);
236  }
237  }
238  hist[nFiles]->SetBinContent(bin, content);
239  hist[nFiles]->SetBinError(bin, error);
240  if (strstr(histName.Data(), "Vtheta")) { // calculate r0theta from Vtheta
241  Vtheta = hist[nFiles]->GetBinContent(bin);
242  VthetaErr = hist[nFiles]->GetBinError(bin);
243  r0V[selN][harN][bin] = Vtheta ? j01 / Vtheta : 0.;
244  r0VErr = Vtheta ? r0V[selN][harN][bin] * VthetaErr / Vtheta : 0.;
245  hist_r0theta[selN][harN]->SetBinContent(bin, r0V[selN][harN][bin]);
246  hist_r0theta[selN][harN]->SetBinError(bin, r0VErr);
247  }
248 
249  if (fromG && strstr(histName.Data(), "ImGtheta")) {
250  reG = histReG[selN][harN][thetaN]->GetBinContent(bin);
251  imG = hist[nFiles]->GetBinContent(bin);
252  Gtheta(reG, imG);
253  histG[selN][harN][thetaN]->SetBinContent(bin, Gtheta.Rho2());
254  }
255 
256  } // bin
257  if (!strstr(histName.Data(), "G")) { break; } // only once if no thetas
258 
259  if (fromG && strstr(histName.Data(), "ReGtheta")) {
260  histReG[selN][harN][thetaN] = (TH1D*)hist[nFiles]->Clone();
261  }
262 
263  if ((fromG) && pageNumber == nNames-1) { // calculate r0theta from Gtheta
264  // Find first minimum of the square of the modulus of G for each theta
265  r0 = hist[nFiles]->GetBinCenter(nBins-1); // default value if no minimum
266  for (int N = 2; N < nBins-2; N++) {
267  G0 = histG[selN][harN][thetaN]->GetBinContent(N);
268  Gnext = histG[selN][harN][thetaN]->GetBinContent(N+1);
269  GnextNext = histG[selN][harN][thetaN]->GetBinContent(N+2);
270 
271  if (Gnext > G0 && GnextNext > G0) { // lowest point
272  Glast = histG[selN][harN][thetaN]->GetBinContent(N-1);
273  Xlast = histG[selN][harN][thetaN]->GetBinCenter(N-1);
274  X0 = histG[selN][harN][thetaN]->GetBinCenter(N);
275  Xnext = histG[selN][harN][thetaN]->GetBinCenter(N+1);
276  r0 = X0;
277  r0 -= ((X0 - Xlast)*(X0 - Xlast)*(G0 - Gnext) - (X0 - Xnext)*(X0 - Xnext)*(G0 - Glast)) /
278  (2.*((X0 - Xlast)*(G0 - Gnext) - (X0 - Xnext)*(G0 - Glast))); // intopolated minimum
279  break;
280  } // if
281  } // N
282  hist_r0theta[selN][harN]->SetBinContent(thetaN+1, r0);
283  hist_r0theta[selN][harN]->SetBinError(thetaN+1, 0.);
284 
285  r0fromV = r0V[selN][harN][thetaN+1];
286 // if (r0fromV != 0.) { cout << "sel, har, theta= " << selN << ", " << harN << ", "
287 // << thetaN << ": " << setprecision(4) << r0fromV << " \t"
288 // << r0 << " \t\t" << setprecision(2)
289 // << 100. * (r0-r0fromV)/r0fromV << "%" << endl; }
290  } // if last page
291  } // thetaN
292  } // harN
293  } // selN
294  } // pageNumber
295 
296  // write the output file
297  //histFile[nFiles]->ls();
298  histFile[nFiles]->Write(0, TObject::kOverwrite);
299  histFile[nFiles]->Close();
300  delete histFile[nFiles];
301 
302  // copy the output back
303  nFile = 0;
304  fstream ascii_out;
305  ascii_out.open(listFileName, ios::in);
306  ascii_out >> rootFileName;
307  while(!ascii_out.eof()) {
308  if (strstr(rootFileName,".root")==0) continue;
309  cout << nFile+1 << " file name = " << rootFileName << endl;
310  sprintf(cpCommand, "cp flow.firstPassLYZ%2d.root %s", cenNo, rootFileName);
311  system(cpCommand);
312  ascii_out >> rootFileName;
313  nFile++;
314  }
315  cout << "####" << endl << endl;
316 
317  // clean up
318  sprintf(rmCommand, "rm %s", listFileName);
319  system(rmCommand);
320  sprintf(rmCommand, "rm flow.firstPassLYZ%2d.root", cenNo);
321  system(rmCommand);
322 
323  } // cenNo
324 }
325 
327 //
328 // $Log: doFlowSumFirstPass.C,v $
329 // Revision 1.2 2006/07/06 16:58:37 posk
330 // Calculation of v1 for LYZ selection=2 is done with mixed harmonics.
331 //
332 // Revision 1.1 2006/03/22 21:57:46 posk
333 // Macro to sum the Generating Functions between the two passes.
334 //
335 //
TDataSet::Clone
virtual TObject * Clone(const char *newname="") const
the custom implementation fo the TObject::Clone
Definition: TDataSet.cxx:308
Generated by   doxygen 1.8.5