StRoot: StFlowAnalysisMaker/minBias.C Source File

00001 
00002 //
00003 // $Id: minBias.C,v 1.15 2007/02/06 19:00:51 posk Exp $
00004 //
00005 // Author:       Art Poskanzer and Alexander Wetzler, Mar 2001
00006 //                 Kirill Filimonov treated the one count case
00007 // Description:  Macro to add histograms together.
00008 //               The v histograms will be added with yield weighting.
00009 //               First file anaXX.root given by first run number XX.
00010 //               Last file anaXX.root given by last run number XX.
00011 //               Output file given by output run number.
00012 //               For 10-50% centrality: minBias(X4,X7,X0)
00013 //
00014 //
00016 #ifndef __CINT__
00017 #include <fstream.h>
00018 #include "TSystem.h"
00019 #include <TFile.h>
00020 #include "TH1.h"
00021 #include "TH2.h"
00022 #include "TProfile.h"
00023 #include "TKey.h"
00024 #include "TObject.h"
00025 #endif
00026 
00027 void minBias(Int_t firstRunNo, Int_t lastRunNo, Int_t outputRunNo=99) {
00028 
00029   const   int nCens = lastRunNo - firstRunNo + 1;
00030   int     nSels = 2;
00031   const   int nHars = 4;
00032   char    fileName[80];
00033   TFile*  histFile[nCens+1];
00034   TH1*    hist[nCens+1];
00035   TH2*    yieldPartHist[nCens];
00036   TH1*    yieldPartHistPt[nCens];
00037   TH1*    yieldPartHistEta[nCens];
00038   Float_t yieldTotal[nCens];
00039   
00040   // names of histograms to be added with weighting
00041   const char* baseName[] = { 
00042     "Flow_Res_",
00043     "Flow_v2D_",   // must be before vEta and vPt
00044     "Flow_vEta_",
00045     "Flow_vPt_",
00046     "Flow_v_",
00047     "FlowLYZ_r0theta_",
00048     "FlowLYZ_Vtheta_",
00049     "FlowLYZ_V_",
00050     "FlowLYZ_vr0_",
00051     "FlowLYZ_vEta_",
00052     "FlowLYZ_vPt_",
00053     "FlowLYZ_v_",
00054     "Flow_v2D_ScalarProd_",
00055     "Flow_vEta_ScalarProd_",
00056     "Flow_vPt_ScalarProd_",
00057     "Flow_v_ScalarProd_",
00058     "Flow_Cumul_vEta_Order2_",
00059     "Flow_Cumul_vPt_Order2_",
00060     "Flow_Cumul_v_Order2_",
00061     "Flow_Cumul_vEta_Order4_",
00062     "Flow_Cumul_vPt_Order4_",
00063     "Flow_Cumul_v_Order4_"
00064   };
00065   const int nNames = sizeof(baseName) / sizeof(char*);
00066 
00067   // open the files
00068   for (int n = 0; n < nCens; n++) {
00069     sprintf(fileName, "ana%2d.root", firstRunNo + n);
00070     cout << " file name = " << fileName << endl;
00071     histFile[n] = new TFile(fileName);
00072     if (!histFile[n]) {
00073       cout << "### Can't find file " << fileName << endl;
00074       return;
00075     }
00076   }
00077   sprintf(fileName, "ana%2d.root", outputRunNo);
00078   cout << " output file name = " << fileName << endl << endl;
00079   histFile[nCens] = new TFile(fileName, "RECREATE");
00080     
00081   // add all histograms with no weighting
00082   TKey*    key;
00083   TObject* obj;
00084   TIter nextkey(histFile[0]->GetListOfKeys());  
00085   while (key = (TKey*)nextkey()) {
00086     histFile[0]->cd();
00087     //key->ls();
00088     obj = key->ReadObj();
00089     const char* objName;
00090     if (obj->InheritsFrom("TH1")) { // TH1 or TProfile
00091       hist[0] = (TH1*)obj;
00092       objName = key->GetName();
00093       cout << "hist name= " << objName << endl;
00094       for (int n = 1; n < nCens; n++) {
00095         hist[1] = (TH1*)histFile[n]->Get(objName);
00096         hist[0]->Add(hist[1]);
00097       }
00098     }
00099     histFile[nCens]->cd();
00100     obj->Write(objName);
00101     delete obj;
00102     obj = NULL;
00103   }
00104    
00105   // get yield histogram
00106   cout<<endl;
00107   for (int n = 0; n < nCens; n++) {
00108     yieldPartHist[n] = dynamic_cast<TH2*>(histFile[n]->Get("Flow_YieldPart2D"));
00109     yieldPartHistPt[n] = dynamic_cast<TH1*>(histFile[n]->Get("FlowLYZ_YieldPartPt"));
00110     yieldPartHistEta[n] = dynamic_cast<TH1*>(histFile[n]->Get("FlowLYZ_YieldPartEta"));
00111     if (yieldPartHistPt[n]) { yieldTotal[n] = yieldPartHistPt[n]->Integral(); }
00112   }
00113 
00114   for (int pageNumber = 0; pageNumber < nNames; pageNumber++ ) {
00115     bool twoD = kFALSE;
00116     if (strstr(baseName[pageNumber],"v2D")) twoD = kTRUE;
00117 
00118     for (int selN = 0; selN < nSels; selN++) {
00119       
00120       // no harmonics
00121       bool noHars = kFALSE;
00122       int nHar = nHars;
00123       if (strstr(baseName[pageNumber],"_v_")!=0 ||
00124           strstr(baseName[pageNumber],"Res")!=0 ||
00125           strstr(baseName[pageNumber],"_V_")!=0 ||
00126           strstr(baseName[pageNumber],"_vr0_")!=0) {
00127         noHars = kTRUE;
00128         nHar = 1;
00129       } 
00130       for (int harN = 0; harN < nHar; harN++) { 
00131         
00132         // construct histName
00133         TString* histName = new TString(baseName[pageNumber]);
00134         *histName += "Sel";
00135         *histName += selN+1;
00136         if (!noHars) {  
00137           *histName += "_Har";
00138           *histName += harN+1;
00139         }
00140                 
00141         // get the histograms
00142         for (int n = 0; n < nCens+1; n++) {
00143           hist[n] = dynamic_cast<TH1*>(histFile[n]->Get(histName->Data()));
00144         }
00145         if (hist[0]) { cout << "hist name= " << histName->Data() << endl; }
00146 
00147         const int lastHist = nCens;
00148         int nBins;      // set by 2D of centrality lastHist
00149         if (!hist[lastHist]) continue;
00150         int xBins = hist[lastHist]->GetNbinsX();
00151         int yBins;
00152         TH1F *yieldY, *yieldPt;
00153         if (twoD) {
00154           yBins = hist[lastHist]->GetNbinsY();
00155           nBins = xBins + (xBins + 2) * yBins;
00156           float yMax  = hist[lastHist]->GetXaxis()->GetXmax();
00157           float yMin  = hist[lastHist]->GetXaxis()->GetXmin();
00158           yieldY = new TH1F("Yield_Y", "Yield_Y", xBins, yMin, yMax);
00159           yieldY->SetXTitle("Rapidity");
00160           yieldY->SetYTitle("Counts");
00161           float ptMax  = hist[lastHist]->GetYaxis()->GetXmax();
00162           yieldPt = new TH1F("Yield_Pt", "Yield_Pt", yBins, 0., ptMax);
00163           yieldPt->SetXTitle("Pt (GeV/c)");
00164           yieldPt->SetYTitle("Counts");
00165         } else {
00166           nBins = xBins + 2;
00167         }
00168         
00169         // loop over the bins
00170         if (strstr(baseName[pageNumber],"Res") ||
00171             strstr(baseName[pageNumber],"Cumul")) { // with error weighting
00172           cout << "  With error weighting" << endl;
00173           float content, error, errorSq, meanContent, meanError, weight;
00174           for (int bin = 0; bin < nBins; bin++) {
00175             meanContent = 0.;
00176             meanError   = 0.;
00177             weight      = 0.;
00178             for (int n = 0; n < nCens; n++) {
00179               if (hist[n]) {
00180                 content = hist[n]->GetBinContent(bin);
00181                 error   = hist[n]->GetBinError(bin);
00182                 errorSq = error * error;
00183                 if (errorSq > 0.) {
00184                   meanContent += content / errorSq;
00185                   weight      += 1. / errorSq;
00186                 }
00187               }
00188             }
00189             if (weight > 0.) {
00190             meanContent /= weight;
00191             meanError = sqrt(1. / weight);
00192             hist[nCens]->SetBinContent(bin, meanContent);
00193             hist[nCens]->SetBinError(bin, meanError);
00194             }
00195           }
00196         } else {                                   // with yield weighting
00197           cout << "  With yield weighting" << endl;
00198           float v, verr, content, error, yield, y, pt;
00199           double vSum, vSum2, error2sum, yieldSum, vRms, verrRms, vSumRms, yieldSumRms, vSumRms2;
00200           for (int bin = 0; bin < nBins; bin++) {
00201             v           = 0.;
00202             verr        = 0.;
00203             vSum        = 0.;
00204             vSum2       = 0.;
00205             content     = 0.;
00206             error       = 0.;
00207             error2sum   = 0.;
00208             yield       = 0.;
00209             yieldSum    = 0.;
00210             vRms        = 0.;
00211             verrRms     = 0.;
00212             vSumRms     = 0.;
00213             yieldSumRms = 0.;
00214             vSumRms2    = 0.;
00215             for (int n = 0; n < nCens; n++) {
00216               if (hist[n]) {
00217                 if (strstr(baseName[pageNumber],"LYZ")) {
00218                   if (strstr(histName->Data(), "vEta")) {
00219                     yield = yieldPartHistEta[n]->GetBinContent(bin);
00220                   } else if (strstr(histName->Data(), "vPt")) {
00221                     yield = yieldPartHistPt[n]->GetBinContent(bin);
00222                   } else {                               // r0theta, Vtheta, _V_, vr0, _v_, Gtheta
00223                     yield = yieldTotal[n];
00224                   }
00225                 } else {        
00226                   if (strstr(histName->Data(),"v2D")) {
00227                     yield = yieldPartHist[n]->GetBinContent(bin);
00228                   } else if (strstr(histName->Data(),"vEta")) {
00229                     yield = yieldPartHist[n]->Integral(bin, bin, 1, yBins);
00230                     if (selN==0 && harN==0) {
00231                       y = yieldPartHist[n]->GetXaxis()->GetBinCenter(bin);
00232                       yieldY->Fill(y, yield);
00233                     }
00234                   } else if (strstr(histName->Data(),"vPt")) {
00235                     yield = yieldPartHist[n]->Integral(1, xBins, bin, bin);
00236                     if (selN==0 && harN==0) {
00237                       pt = yieldPartHist[n]->GetYaxis()->GetBinCenter(bin);
00238                       yieldPt->Fill(pt, yield);
00239                     }
00240                   } else {                                        // _v
00241                     yield = yieldPartHist[n]->Integral();
00242                   }
00243                 }
00244                 v = hist[n]->GetBinContent(bin);
00245                 if (yield==1) { // special case to calculate the correct error
00246                   vSumRms     += v;
00247                   yieldSumRms += yield;
00248                   vSumRms2    += v*v;
00249                 } else {
00250                   verr = hist[n]->GetBinError(bin);
00251                   if (v != 0 ) {
00252                     yieldSum  += yield;
00253                     vSum      += yield * v;
00254                     vSum2     += v * v * yield;
00255                     error2sum += pow(yield * verr, 2.);
00256                   }
00257                 }
00258               }
00259             }
00260             
00261             if (yieldSumRms) vRms = vSumRms / yieldSumRms;
00262             if (yieldSumRms>1) {
00263               verrRms    = sqrt(vSumRms2 - vSumRms*vSumRms / yieldSumRms)
00264                 / yieldSumRms;
00265               yieldSum  += yieldSumRms;
00266               vSum      += vSumRms;
00267               error2sum += pow(yieldSumRms * verrRms, 2.);
00268               vSum2     += vRms * vRms * yieldSumRms;
00269             }
00270             
00271             if (yieldSum) {
00272               content = vSum / yieldSum;
00273               if (yieldSumRms==1) {
00274                 error = sqrt(error2sum + vSum2 - vSum*vSum / yieldSum) / yieldSum;
00275                 error = yieldSum / (yieldSum+1) * sqrt(error*error + 
00276                    (vRms - content)*(vRms - content) / (yieldSum * (yieldSum+1)));
00277               } else {
00278                 error = sqrt(error2sum + vSum2 - vSum*vSum/yieldSum) / yieldSum;
00279               }
00280               hist[nCens]->SetBinContent(bin, content);
00281               hist[nCens]->SetBinError(bin, error);
00282             }
00283           }
00284         } 
00285         delete histName;
00286       }
00287     }
00288   }
00289   
00290   //histFile[nCens]->ls();
00291   histFile[nCens]->Write(0, TObject::kOverwrite);
00292   histFile[nCens]->Close();
00293   delete histFile[nCens];
00294   
00295 }
00296 
00298 //
00299 // $Log: minBias.C,v $
00300 // Revision 1.15  2007/02/06 19:00:51  posk
00301 // In Lee Yang Zeros method, introduced recentering of Q vector.
00302 // Reactivated eta symmetry cut.
00303 //
00304 // Revision 1.13  2006/03/22 22:02:07  posk
00305 // Updates to macros.
00306 //
00307 // Revision 1.12  2006/02/24 18:36:04  posk
00308 // Updated for LeeYangZerosMaker histograms.
00309 //
00310 // Revision 1.11  2004/03/01 22:43:42  posk
00311 // Changed some "->" to ".".
00312 //
00313 // Revision 1.10  2003/08/26 21:10:12  posk
00314 // Calculates v8 if nHars=8.
00315 //
00316 // Revision 1.9  2003/06/27 21:25:44  posk
00317 // v4 and v6 are with repect to the 2nd harmonic event plane.
00318 //
00319 // Revision 1.8  2003/03/11 23:03:07  posk
00320 // Includes scalar product and cumulant hists.
00321 //
00322 // Revision 1.7  2002/06/11 21:54:15  posk
00323 // Kirill's further correction to minBias.C for bins with one count.
00324 //
00325 // Revision 1.6  2002/05/21 18:42:18  posk
00326 // Kirill's correction to minBias.C for bins with one count.
00327 //
00328 // Revision 1.5  2001/11/09 21:15:00  posk
00329 // Switched from CERNLIB to TMath. Using global dca instead of dca.
00330 //
00331 // Revision 1.4  2001/05/22 20:05:47  posk
00332 // Now outputs a hist.root file.
00333 // The v values are averaged with yield weighting.
00334 //
00335 // Revision 1.3  2000/09/29 22:53:17  posk
00336 // More histograms.
00337 //
00338 // Revision 1.2  2000/09/26 20:54:11  posk
00339 // Updated documentation.
00340 //
00341 // Revision 1.1  2000/09/26 00:19:43  posk
00342 // New macro to add centrality-selected histograms with proper weights, to make
00343 // minimum bias histogram.
00344 //
00345 //