StRoot: StFlowAnalysisMaker/doFlowSumFirstPass.C Source File

00001 
00002 //
00003 // $Id: doFlowSumFirstPass.C,v 1.2 2006/07/06 16:58:37 posk Exp $
00004 //
00005 // Makes temporary root.files.<cenNo> containing lists of flow.firstPassLYZNew.root files
00006 //   in all subdirectories of outDir.
00007 // outDir is a link in this directory.
00008 // The cenNo is the firstCenNo + cen.
00009 // The link argument is the prefix of the subdirectory name.
00010 // Adds the histograms together by yield weighted averaging for all files for this centrality.
00011 // Writes out files called flow.firstPassLYZ<cenNo>.root.
00012 // Copies them back to the subdirectories as flow.firstPassLYZNew.root.
00013 // This macro must be in your directory containing the outDir link
00014 //   and run between the first and second passes.
00015 //
00016 // by Art Poskanzer, Feb 2006
00017 //
00019 
00020 #include <iostream.h>
00021 #include <iomanip.h>
00022 #include <fstream.h>
00023 #include "TH1.h"
00024 #include "TFile.h"
00025 #include "TComplex.h"
00026 
00027 void doFlowSumFirstPass(char*  link, int firstCenNo) {
00028 
00029   int    nCens = 9; // 9
00030   Bool_t fromG = kTRUE; // the prefered setting
00031 
00032   char    rootFileName[80];
00033   char    rootDirName[80];
00034   char    listFileName[80];
00035   char    rootOutName[80];
00036   char    lsCommand[80];
00037   char    firstPassCommand[80];
00038   char    rmCommand[80];
00039   char    cpCommand[80];
00040   TFile*  histFile[1000];
00041   Float_t j01 = 2.405;
00042   int     maxTheta;
00043   const   int nSels = 2;
00044   const   int nHars = 2;
00045 
00046   // names of histograms to be averaged with weighting (no profiles and no 2D hists)
00047   const char* baseName[] = { 
00048     "FlowLYZ_r0theta_",
00049     "FlowLYZ_Vtheta_", // must come after r0theta
00050     "FlowLYZ_Gtheta",  // first G
00051     "FlowReGtheta",
00052     "FlowImGtheta"     // after ReG
00053   };
00054   const int nNames = sizeof(baseName) / sizeof(char*);
00055   
00056   for (int cen = 0; cen < nCens; cen++) {
00057     int cenNo = firstCenNo + cen;
00058 
00059     // open the input files
00060     Int_t nFile = 0;
00061     sprintf(listFileName, "root.files.%2d", cenNo);
00062     sprintf(lsCommand, "ls -1 outDir/%s-*-%2d/flow.firstPassLYZNew.root > %s", link, cenNo,
00063             listFileName);
00064     //sprintf(lsCommand, "ls -1 flow.firstPassLYZNew.root > %s", listFileName);
00065     system(lsCommand);
00066     fstream ascii_in;
00067     ascii_in.open(listFileName, ios::in);
00068     ascii_in >> rootFileName;
00069     while(!ascii_in.eof()) {
00070       if (strstr(rootFileName,".root")==0) continue;
00071       histFile[nFile] = new TFile(rootFileName);
00072       //if (nFile==0) { histFile[nFile]->ls(); }
00073       strcpy(rootDirName, rootFileName);
00074       char* cp = strstr(rootDirName, "flow.");
00075       if (cp) *cp = '\0';                    // truncate to directory name
00076       sprintf(firstPassCommand, "test -f %s", rootFileName);
00077       if (system(firstPassCommand)) {
00078         cout << "### No first pass for " << rootDirName << endl;
00079         //delete histFile[nFile];
00080         continue;
00081       }
00082       cout << nFile+1 << " directory name = " << rootDirName << endl;
00083       ascii_in >> rootFileName;
00084       nFile++;
00085     }
00086     const Int_t nFiles = nFile;
00087     cout << "input files: " << nFiles << endl << endl;
00088     if (nFiles == 0) { return; }
00089     
00090     // the following 3 lines will prevent compiling this macro
00091     TH1D*    hist[nFiles+1];  // last file is the output
00092     TH1D*    hist_r0theta[nSels][nHars];
00093     TH1F*    yieldMultHist[nFiles+1];
00094     Float_t  yieldTotal[nFiles];
00095     
00096     // open the output file
00097     sprintf(rootOutName, "flow.firstPassLYZ%2d.root", cenNo);
00098     cout << "tempory output file name = " << rootOutName << endl << endl ;
00099     histFile[nFiles] = new TFile(rootOutName, "RECREATE"); // last file is new out file
00100         
00101     // get multiplicity histograms
00102     for (int n = 0; n < nFiles; n++) {
00103       yieldMultHist[n] = dynamic_cast<TH1F*>(histFile[n]->Get("FlowLYZ_Mult"));
00104       if (!yieldMultHist[n]) {
00105         cout << "### dynamic cast can't find histogram FlowLYZ_Mult" << endl;
00106         return;
00107       }
00108       yieldTotal[n] = yieldMultHist[n]->GetEntries() * yieldMultHist[n]->GetMean();
00109       cout << n+1 << ": " << setprecision(7) << yieldTotal[n] << " particles" << endl;
00110     }
00111 
00112     // add mult hists together for output
00113     histFile[nFiles]->cd();
00114     yieldMultHist[nFiles] = (TH1F*)yieldMultHist[0]->Clone();
00115     for (int n = 1; n < nFiles; n++) {
00116       yieldMultHist[nFiles]->Add(yieldMultHist[n]);
00117     }
00118     cout << endl;
00119 
00120     // Get maxTheta
00121     if (cen==0) {
00122       TString histNameNtheta("FlowLYZ_r0theta_Sel1_Har2");
00123       TH1D* histNtheta = (TH1D*)histFile[0]->Get(histNameNtheta);
00124       if (!histNtheta) {
00125         cout << "### Can't find file " << histNameNtheta << endl;
00126         return;
00127       }
00128       maxTheta = histNtheta->GetNbinsX();
00129       cout << "maxTheta = " << maxTheta << endl << endl;
00130       const int thetas = maxTheta + 2;
00131     }
00132     TH1D*    histReG[nSels][nHars][thetas];
00133     TH1D*    histG[nSels][nHars][thetas];
00134 
00135     //if (fromG) { cout << "r0:\t\t\t from V \t from G \t change" << endl; }
00136     float r0V[nSels][nHars][thetas];
00137     float Xlast, X0, Xnext;
00138     double reG, imG, Glast, G0, Gnext, GnextNext;
00139     TComplex Gtheta;
00140 
00141     // with yield weighted averaging
00142     for (int pageNumber = 0; pageNumber < nNames; pageNumber++ ) {
00143       for (int selN = 0; selN < nSels; selN++) { 
00144         for (int harN = 0; harN < nHars; harN++) {
00145           for (int thetaN = 0; thetaN < maxTheta; thetaN++) {   
00146           
00147             // construct histName
00148             TString histName(baseName[pageNumber]);
00149             if (strstr(histName.Data(), "G")) {
00150               histName += thetaN;
00151               histName += "_Sel";
00152               histName += selN+1;
00153               histName += "_Har";
00154               histName += harN+1;
00155             } else {
00156               histName += "Sel";
00157               histName += selN+1;
00158               histName += "_Har";
00159               histName += harN+1;
00160             }
00161             //cout << "hist name= " << histName << endl;
00162             
00163             // get the histograms
00164             for (int n = 0; n < nFiles; n++) {
00165               hist[n] = dynamic_cast<TH1D*>(histFile[n]->Get(histName));
00166               if (!hist[n]) {
00167                 cout << "### dynamic cast can't find file " << n << ": histogram " <<
00168                   histName << endl;
00169                 return;
00170               }
00171             }
00172             
00173             // make the output hists
00174             histFile[nFiles]->cd();
00175             hist[nFiles] = new TH1D(*(hist[0]));
00176             
00177             // get pointer to the r0 output hist
00178             if (strstr(histName.Data(), "r0theta")) {
00179               hist_r0theta[selN][harN] = dynamic_cast<TH1D*>(histFile[nFiles]->Get(histName));
00180               if (!hist_r0theta[selN][harN]) {
00181                 cout << "### dynamic castan't find hist " << histName << endl;
00182                 return;
00183               }
00184             }
00185             
00186             // get pointer to the G output hist
00187             if (strstr(histName.Data(), "_G")) {
00188               histG[selN][harN][thetaN] = dynamic_cast<TH1D*>(histFile[nFiles]->Get(histName));
00189               if (!histG[selN][harN][thetaN]) {
00190                 cout << "### dynamic cast can't find hist " << histName << endl;
00191                 return;
00192               }
00193             }
00194 
00195             // get the number of bins for this hist
00196             if (!hist[0]) continue;
00197             int nBins = hist[0]->GetNbinsX() + 2;
00198             
00199             // loop over the bins
00200             float  yield, r0fromV;
00201             double content, error, Vtheta, VthetaErr, r0, r0VErr;
00202             double v, vSum, v2Sum, error2Sum, yieldSum, yieldSum2, yield2Sum;
00203             for (int bin = 0; bin < nBins; bin++) {
00204               v         = 0.;
00205               vSum      = 0.;
00206               v2Sum     = 0.;
00207               content   = 0.;
00208               error     = 0.;
00209               error2Sum = 0.;
00210               yield     = 0.;
00211               yieldSum  = 0.;
00212               yield2Sum = 0.;
00213               for (int n = 0; n < nFiles; n++) { // loop over the files
00214                 if (hist[n]) {
00215                   yield = yieldTotal[n];
00216                   v = hist[n]->GetBinContent(bin);
00217                   if (v != 0. && yield > 1.) { // error is wrong for one count
00218                     yieldSum  += yield;
00219                     yield2Sum += yield*yield;
00220                     vSum      += yield * v;
00221                     v2Sum     += yield*yield * v*v;
00222                     error2Sum += pow(yield * hist[n]->GetBinError(bin), 2.);
00223                   }
00224                 }
00225               } // nFiles
00226               if (yieldSum) {
00227                 content   = vSum / yieldSum;
00228                 yieldSum2 = yieldSum*yieldSum;
00229                 v2Sum     /= yieldSum2;
00230                 yield2Sum /= yieldSum2;
00231                 error2Sum /= yieldSum2;
00232                 if (strstr(histName.Data(), "G")) {
00233                   error = 0.;
00234                 } else { 
00235                   error = sqrt(error2Sum);
00236                 }
00237               }
00238               hist[nFiles]->SetBinContent(bin, content);
00239               hist[nFiles]->SetBinError(bin, error);
00240               if (strstr(histName.Data(), "Vtheta")) { // calculate r0theta from Vtheta
00241                 Vtheta = hist[nFiles]->GetBinContent(bin);
00242                 VthetaErr = hist[nFiles]->GetBinError(bin);
00243                 r0V[selN][harN][bin] = Vtheta ? j01 / Vtheta : 0.;
00244                 r0VErr = Vtheta ? r0V[selN][harN][bin] * VthetaErr / Vtheta : 0.;
00245                 hist_r0theta[selN][harN]->SetBinContent(bin, r0V[selN][harN][bin]);
00246                 hist_r0theta[selN][harN]->SetBinError(bin, r0VErr);
00247               }
00248 
00249               if (fromG && strstr(histName.Data(), "ImGtheta")) {
00250                 reG = histReG[selN][harN][thetaN]->GetBinContent(bin);
00251                 imG = hist[nFiles]->GetBinContent(bin);
00252                 Gtheta(reG, imG);
00253                 histG[selN][harN][thetaN]->SetBinContent(bin, Gtheta.Rho2());
00254               }
00255 
00256             } // bin
00257             if (!strstr(histName.Data(), "G")) { break; } // only once if no thetas
00258 
00259             if (fromG && strstr(histName.Data(), "ReGtheta")) {
00260               histReG[selN][harN][thetaN] = (TH1D*)hist[nFiles]->Clone();
00261             }
00262 
00263             if ((fromG) && pageNumber == nNames-1) { // calculate r0theta from Gtheta
00264               // Find first minimum of the square of the modulus of G for each theta
00265               r0 = hist[nFiles]->GetBinCenter(nBins-1); // default value if no minimum
00266               for (int N = 2; N < nBins-2; N++) {
00267                 G0        = histG[selN][harN][thetaN]->GetBinContent(N);
00268                 Gnext     = histG[selN][harN][thetaN]->GetBinContent(N+1);
00269                 GnextNext = histG[selN][harN][thetaN]->GetBinContent(N+2);
00270                 
00271                 if (Gnext > G0 && GnextNext > G0) { // lowest point
00272                   Glast     = histG[selN][harN][thetaN]->GetBinContent(N-1);
00273                   Xlast     = histG[selN][harN][thetaN]->GetBinCenter(N-1);
00274                   X0        = histG[selN][harN][thetaN]->GetBinCenter(N);
00275                   Xnext     = histG[selN][harN][thetaN]->GetBinCenter(N+1);
00276                   r0 = X0;
00277                   r0 -= ((X0 - Xlast)*(X0 - Xlast)*(G0 - Gnext) - (X0 - Xnext)*(X0 - Xnext)*(G0 - Glast)) /
00278                     (2.*((X0 - Xlast)*(G0 - Gnext) - (X0 - Xnext)*(G0 - Glast))); // intopolated minimum
00279                   break;
00280                 } // if
00281               } // N
00282               hist_r0theta[selN][harN]->SetBinContent(thetaN+1, r0);
00283               hist_r0theta[selN][harN]->SetBinError(thetaN+1, 0.);
00284 
00285               r0fromV = r0V[selN][harN][thetaN+1];
00286 //            if (r0fromV != 0.) { cout << "sel, har, theta= " << selN << ", " << harN << ", "
00287 //                                      << thetaN << ": " << setprecision(4) << r0fromV << " \t"
00288 //                                      << r0 << " \t\t" << setprecision(2)
00289 //                                      << 100. * (r0-r0fromV)/r0fromV << "%" << endl; }
00290             } // if last page
00291           } // thetaN
00292         } // harN
00293       } // selN
00294     } // pageNumber
00295 
00296     // write the output file
00297     //histFile[nFiles]->ls();
00298     histFile[nFiles]->Write(0, TObject::kOverwrite);
00299     histFile[nFiles]->Close();
00300     delete histFile[nFiles];
00301 
00302     // copy the output back
00303     nFile = 0;
00304     fstream ascii_out;
00305     ascii_out.open(listFileName, ios::in);
00306     ascii_out >> rootFileName;
00307     while(!ascii_out.eof()) {
00308       if (strstr(rootFileName,".root")==0) continue;
00309       cout << nFile+1 << " file name = " << rootFileName << endl;
00310       sprintf(cpCommand, "cp flow.firstPassLYZ%2d.root %s", cenNo, rootFileName);
00311       system(cpCommand);
00312       ascii_out >> rootFileName;
00313       nFile++;
00314     }
00315     cout << "####" << endl << endl;
00316 
00317     // clean up
00318     sprintf(rmCommand, "rm %s", listFileName);
00319     system(rmCommand);
00320     sprintf(rmCommand, "rm flow.firstPassLYZ%2d.root", cenNo);
00321     system(rmCommand);
00322     
00323   } // cenNo
00324 }
00325 
00327 //
00328 // $Log: doFlowSumFirstPass.C,v $
00329 // Revision 1.2  2006/07/06 16:58:37  posk
00330 // Calculation of v1 for LYZ selection=2 is done with mixed harmonics.
00331 //
00332 // Revision 1.1  2006/03/22 21:57:46  posk
00333 // Macro to sum the Generating Functions between the two passes.
00334 //
00335 //