27 void doFlowSumFirstPass(
char* link,
int firstCenNo) {
32 char rootFileName[80];
34 char listFileName[80];
37 char firstPassCommand[80];
40 TFile* histFile[1000];
47 const char* baseName[] = {
54 const int nNames =
sizeof(baseName) /
sizeof(
char*);
56 for (
int cen = 0; cen < nCens; cen++) {
57 int cenNo = firstCenNo + cen;
61 sprintf(listFileName,
"root.files.%2d", cenNo);
62 sprintf(lsCommand,
"ls -1 outDir/%s-*-%2d/flow.firstPassLYZNew.root > %s", link, cenNo,
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);
73 strcpy(rootDirName, rootFileName);
74 char* cp = strstr(rootDirName,
"flow.");
76 sprintf(firstPassCommand,
"test -f %s", rootFileName);
77 if (system(firstPassCommand)) {
78 cout <<
"### No first pass for " << rootDirName << endl;
82 cout << nFile+1 <<
" directory name = " << rootDirName << endl;
83 ascii_in >> rootFileName;
86 const Int_t nFiles = nFile;
87 cout <<
"input files: " << nFiles << endl << endl;
88 if (nFiles == 0) {
return; }
92 TH1D* hist_r0theta[nSels][nHars];
93 TH1F* yieldMultHist[nFiles+1];
94 Float_t yieldTotal[nFiles];
97 sprintf(rootOutName,
"flow.firstPassLYZ%2d.root", cenNo);
98 cout <<
"tempory output file name = " << rootOutName << endl << endl ;
99 histFile[nFiles] =
new TFile(rootOutName,
"RECREATE");
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;
108 yieldTotal[n] = yieldMultHist[n]->GetEntries() * yieldMultHist[n]->GetMean();
109 cout << n+1 <<
": " << setprecision(7) << yieldTotal[n] <<
" particles" << endl;
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]);
122 TString histNameNtheta(
"FlowLYZ_r0theta_Sel1_Har2");
123 TH1D* histNtheta = (TH1D*)histFile[0]->Get(histNameNtheta);
125 cout <<
"### Can't find file " << histNameNtheta << endl;
128 maxTheta = histNtheta->GetNbinsX();
129 cout <<
"maxTheta = " << maxTheta << endl << endl;
130 const int thetas = maxTheta + 2;
132 TH1D* histReG[nSels][nHars][thetas];
133 TH1D* histG[nSels][nHars][thetas];
136 float r0V[nSels][nHars][thetas];
137 float Xlast, X0, Xnext;
138 double reG, imG, Glast, G0, Gnext, GnextNext;
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++) {
148 TString histName(baseName[pageNumber]);
149 if (strstr(histName.Data(),
"G")) {
164 for (
int n = 0; n < nFiles; n++) {
165 hist[n] =
dynamic_cast<TH1D*
>(histFile[n]->Get(histName));
167 cout <<
"### dynamic cast can't find file " << n <<
": histogram " <<
174 histFile[nFiles]->cd();
175 hist[nFiles] =
new TH1D(*(hist[0]));
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;
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;
196 if (!hist[0])
continue;
197 int nBins = hist[0]->GetNbinsX() + 2;
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++) {
213 for (
int n = 0; n < nFiles; n++) {
215 yield = yieldTotal[n];
216 v = hist[n]->GetBinContent(bin);
217 if (v != 0. && yield > 1.) {
219 yield2Sum += yield*yield;
221 v2Sum += yield*yield * v*v;
222 error2Sum += pow(yield * hist[n]->GetBinError(bin), 2.);
227 content = vSum / yieldSum;
228 yieldSum2 = yieldSum*yieldSum;
230 yield2Sum /= yieldSum2;
231 error2Sum /= yieldSum2;
232 if (strstr(histName.Data(),
"G")) {
235 error = sqrt(error2Sum);
238 hist[nFiles]->SetBinContent(bin, content);
239 hist[nFiles]->SetBinError(bin, error);
240 if (strstr(histName.Data(),
"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);
249 if (fromG && strstr(histName.Data(),
"ImGtheta")) {
250 reG = histReG[selN][harN][thetaN]->GetBinContent(bin);
251 imG = hist[nFiles]->GetBinContent(bin);
253 histG[selN][harN][thetaN]->SetBinContent(bin, Gtheta.Rho2());
257 if (!strstr(histName.Data(),
"G")) {
break; }
259 if (fromG && strstr(histName.Data(),
"ReGtheta")) {
260 histReG[selN][harN][thetaN] = (TH1D*)hist[nFiles]->
Clone();
263 if ((fromG) && pageNumber == nNames-1) {
265 r0 = hist[nFiles]->GetBinCenter(nBins-1);
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);
271 if (Gnext > G0 && GnextNext > G0) {
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);
277 r0 -= ((X0 - Xlast)*(X0 - Xlast)*(G0 - Gnext) - (X0 - Xnext)*(X0 - Xnext)*(G0 - Glast)) /
278 (2.*((X0 - Xlast)*(G0 - Gnext) - (X0 - Xnext)*(G0 - Glast)));
282 hist_r0theta[selN][harN]->SetBinContent(thetaN+1, r0);
283 hist_r0theta[selN][harN]->SetBinError(thetaN+1, 0.);
285 r0fromV = r0V[selN][harN][thetaN+1];
298 histFile[nFiles]->Write(0, TObject::kOverwrite);
299 histFile[nFiles]->Close();
300 delete histFile[nFiles];
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);
312 ascii_out >> rootFileName;
315 cout <<
"####" << endl << endl;
318 sprintf(rmCommand,
"rm %s", listFileName);
320 sprintf(rmCommand,
"rm flow.firstPassLYZ%2d.root", cenNo);
virtual TObject * Clone(const char *newname="") const
the custom implementation fo the TObject::Clone