PWGCF/Correlations/DPhi/FourierDecomposition/lrc/SumAndProjectTH2s.C

   1 // Useful on the prompt:
   2 // TH2* h = (TH2*)Cont_cloizide_DhcAna->FindObject("hS482"); h->Draw("surf1"); h->Integral()
   3 // To quickly count events:
   4 // TH2* h = (TH2*)Cont_cloizide_DhcAna->FindObject("fHEvt");
   5 // h->Draw("colz"); h->Integral()
   6 // TH1D* hc = h->ProjectionY("hc"); hc->Draw();
   7 // hc->Integral(1, 90)
   8
   9 bool ispp = 1;
  10
  11 #include "TH1.h"
  12 #include "TH2.h"
  13 #include "TF1.h"
  14 #include "TFile.h"
  15 #include "TGraphErrors.h"
  16 #include "TCanvas.h"
  17 #include "TString.h"
  18 #include "TLatex.h"
  19 #include "TSystem.h"
  20 #include "TTimeStamp.h"
  21 #include "TStopwatch.h"
  22 #include "TStyle.h"
  23 #include "TMath.h"
  24 #include "TKey.h"
  25 #include "TROOT.h"
  26 #include "TCutG.h"
  27 #include "TFormula.h"
  28 #include "TProfile2D.h"
  29
  30 TFile* inFile  = new TFile("$MYDATA/lhc11a_28nov2011_sum.root", "read");
  31 TFile* outFile = 0;
  32 TString listName = "Cont_cloizide_DhcAna";
  33
  34 // histograms
  35 const int maxHists  = 99999;
  36 TH2 *sHists[maxHists];
  37 TH2 *mHists[maxHists];
  38 TH2 *cHists[maxHists];
  39 TH2 *hEvt  = 0;
  40 TH2 *hTrk  = 0;
  41 TH1 *hDefi = 0; //fHPtTrg
  42 TH1 *hDefj = 0; //fHPtAss
  43 TH1 *hDefk = 0; //fHCent
  44 TH1 *hDefz = 0; //fHZvtx
  45 TFormula *fIndex = 0;
  46
  47 // The angular regions defining the 1d projections. The arrays
  48 // k{Phi,Eta}{Min,Max}[] define these regions.
  49 const int kNRegions = 5;
  50 const char* regionStr[] =
  51   {"NSJET", "RIDGE", "ALL", "ETA_NS", "ETA_AS"};
  52 enum eRegion {NSJET, RIDGE, ALL, ETA_NS, ETA_AS};
  53
  54 // Correlation types:
  55 // 0. Same
  56 // 1. Mixed
  57 // 2. CF def "A" (proj 2Ds to dphi then divide)
  58 // 3. CF def "B" (divide 2D's, then project to dhpi). The def. B CFs
  59 // already exist in the input file, just pass them through to the new
  60 // output file, renaming "c" --> "cB".
  61 const int nCorrTypes = 4;
  62 enum eCorrType {kS=0, kM=1, kCA=2, kCB=3};
  63 TString sCorrType[] = {"s","m","c","c"};
  64
  65 // There are two sets of centrality binning: the cbh bins are the
  66 // original ones from the task. The cb1,2 bins are a superset
  67 // including cbh1,2 + any desired combinations. In Loop 2 (the
  68 // cent. combination loop), cb1,2 are indexed as "k", while cbh1,2 are
  69 // indexed as "cb". If no combined bins are needed, Loop 2 just copies
  70 // the existing bins to the output file.
  71
  72 // Binning for gsi03+ generation (TH2 train output)
  73 int nCentBins = 19; // 10 orig + 8 new combined bins
  74 double cb1[] =  {0, 0, 0,2,2, 1,3,0,    0, 1, 2, 3, 4,  5, 10, 20, 30, 40, 60 };
  75 double cb2[] =  {10,20,2,5,10,3,5,5,    1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 90 };
  76
  77 // Just for reference, this should be equivalent to cbh1,2[] in gsi03:
  78 // double cb1[] =  { 0, 1, 2, 3, 4,  5, 10, 20, 30, 40 };
  79 // double cb2[] =  { 1, 2, 3, 4, 5, 10, 20, 30, 40, 50 };
  80
  81 int nZbins = 0;        // zvtx binning from task - read from input file
  82 double z1[99]   = { 0 };
  83 double z2[99]   = { 0 };
  84 int nCentBinsHist = 0; // The exact binning from task - read from input file
  85 double cbh1[99] =  { 0 };
  86 double cbh2[99] =  { 0 };
  87
  88 int colors[] = {kYellow+2, kOrange+2, kBlack, kMagenta, kRed, kOrange,
  89                 kGreen+2, kCyan, kBlue, kViolet, kGray};
  90
  91 // TGraphs to contain complete bin information. Convention: i,j,k,z
  92 // always means pt_trig, pt_assoc, centrality, and z_vtx bin index,
  93 // respectively.
  94 TGraphErrors *gi, *gj, *gk, *gz;
  95 TObjArray* hists = new TObjArray();
  96 TObjArray* profs = new TObjArray(); // profiles storing <pt> info
  97
  98 // Histos used to compute weight factors for summation.
  99 TH1D* hz[101]; // z-vertex dist. at each cent bin
 100 TH1D* hcent = 0; // No longer used - can remove.
 101
 102 double kPhiMin[kNRegions];
 103 double kPhiMax[kNRegions];
 104 double kEtaMin[kNRegions];
 105 double kEtaMax[kNRegions];
 106
 107 // Reassign as fn. arguments
 108 Double_t dEtaOuter = 1.799;
 109 Double_t dEtaInner = 0.8;
 110
 111 // Functions ---------------------------------------------------------
 112 const char* Region(int r);
 113 char* BinLabel(int ijk, int i);
 114 TGraphErrors* BinGraph(int ijk);
 115 void PrintBins();
 116 TObjArray* GetHistList(TFile& fromfile, TString clname, TString dir);
 117 void MakeWeightHistos();
 118 TH1* ProjectTH2(TH2& h, TString xy, float x1, float y1, float x2, float y2,
 119            TString name, TString opt);
 120 void LoadHistos();
 121 TH2 *GetHist(int i, int j, int z, int k, char type);
 122
 123 // Implementation ----------------------------------------------------
 124 void MakeWeightHistos()
 125 {
 126   // make a z-vertex distribution for each centrality bin
 127   for (int cb=0; cb<nCentBinsHist; cb++) {
 128     TAxis* ax = hEvt->GetYaxis();
 129     int bin1 = ax->FindBin(cbh1[cb]);
 130     int bin2 = ax->FindBin(cbh2[cb]-0.001);
 131
 132     hz[cb] = hEvt->ProjectionX(Form("hz%d", cb), bin1, bin2);
 133
 134     if (0)
 135       hz[cb]->Draw(cb==0? "" : "same");
 136
 137     // now filled from file
 138     // hcent->SetBinContent(cb+1, hz[cb]->Integral());
 139     // hcent->GetXaxis()->SetBinLabel(cb+1,Form("%.2g-%.2g", cbh1[cb],cbh2[cb]));
 140   }
 141
 142   //gk->Draw("aep");
 143   return;
 144 }
 145
 146 void LoadHistos()
 147 {
 148   if (!inFile) {
 149     Error("LoadHistos()","TFile ptr zero");
 150     gSystem->Exit(123);
 151   }
 152
 153   TList *list = 0;
 154   if (ispp)
 155     list = (TList*)inFile->Get("Cont_loizides_DhcAna");
 156   else
 157     list = (TList*)inFile->Get("Cont_cloizide_DhcAna");
 158   if (!list) {
 159     Error("LoadHistos()","TList ptr zero");
 160     gSystem->Exit(123);
 161   }
 162
 163   hEvt = (TH2*)list->FindObject("fHEvt");
 164   hTrk = (TH2*)list->FindObject("fHTrk");
 165   if (!hEvt)
 166     Error("LoadHistos()", "!hEvt");
 167   if (!hTrk)
 168     Error("LoadHistos()", "!hTrk");
 169
 170   hDefi   = (TH1*)list->FindObject("fHPtTrg");
 171   hDefj   = (TH1*)list->FindObject("fHPtAss");
 172   hDefk   = (TH1*)list->FindObject("fHCent");
 173   hDefz   = (TH1*)list->FindObject("fHZvtx");
 174
 175   gi = new TGraphErrors();
 176   gj = new TGraphErrors();
 177   gk = new TGraphErrors();
 178   gz = new TGraphErrors();
 179   gi->SetName("TrigPtBins");
 180   gj->SetName("AsscPtBins");
 181   gk->SetName("EvCentBins");
 182   gz->SetName("EvZvtxBins");
 183
 184   for (int i=0; i<hDefi->GetNbinsX(); i++) {
 185     gi->SetPoint(i, hDefi->GetBinCenter(i+1), 1.0);
 186     gi->SetPointError(i, 0.5*hDefi->GetBinWidth(i+1), 0);
 187   }
 188   for (int i=0; i<hDefj->GetNbinsX(); i++) {
 189     gj->SetPoint(i, hDefj->GetBinCenter(i+1), 1.0);
 190     gj->SetPointError(i, 0.5*hDefj->GetBinWidth(i+1), 0);
 191   }
 192   nCentBinsHist = hDefk->GetNbinsX();
 193   for (int i=0; i<hDefk->GetNbinsX(); i++) {
 194     cbh1[i] = hDefk->GetBinLowEdge(i+1);
 195     cbh2[i] = hDefk->GetBinLowEdge(i+2);
 196   }
 197
 198   // Get TProfiles with <pt> and <pt^2> info into memory
 199   for (int n=0; n<list->GetEntries(); n++) {
 200     TObject* obj = list->At(n);
 201     TString clName = obj->ClassName();
 202     if (clName.Contains("TProfile")) {
 203       profs->Add(obj);
 204     }
 205   }
 206
 207   //  hcent = new TH1D("hcent", "hcent", nCentBinsHist,0,nCentBinsHist);
 208
 209   for (int i=0; i<nCentBins; i++) {
 210     gk->SetPoint(i, (cb1[i]+cb2[i])/2., 1.0);
 211     gk->SetPointError(i, TMath::Abs(cb2[i]-cb1[i])/2, 0);
 212     //cout << i << " " <<  (cb1[i]+cb2[i])/2. << endl;
 213   }
 214
 215   nZbins = hDefz->GetNbinsX();
 216   for (int i=0; i<hDefz->GetNbinsX(); i++) {
 217     gz->SetPoint(i, hDefz->GetBinCenter(i+1), 1.0);
 218     gz->SetPointError(i, 0.5*hDefz->GetBinWidth(i+1), 0);
 219     z1[i] = hDefz->GetBinLowEdge(i+1);
 220     z2[i] = hDefz->GetBinLowEdge(i+2);
 221   }
 222
 223   fIndex = new TFormula("GlobIndex", //this is now relative to 0 (and not to 1 as for histos)
 224                         "(t)*[0]*[1]*[2]+(z)*[0]*[1]+(x)*[0]+(y)+0*[4]");
 225   fIndex->SetParameters(gi->GetN(),
 226                         gj->GetN(),
 227                         gz->GetN(),
 228                         gk->GetN());
 229   fIndex->SetParNames("NTrigBins (i)","NAssocBins (j)", "NZvertexBins (z)", "NCentBins (k)");
 230
 231   for (int i=0;i<maxHists;++i) {
 232     mHists[i]=0;
 233     sHists[i]=0;
 234   }
 235
 236   Int_t ent=list->GetEntries();
 237   for (Int_t i=0;i<ent;++i) {
 238     TH2 *obj = dynamic_cast<TH2*>(list->At(i));
 239     if (!obj)
 240       continue;
 241     TString name(obj->GetName());
 242     if (!name.BeginsWith("hM") && !name.BeginsWith("hS"))
 243       continue;
 244     const char *ptr = name.Data()+2;
 245     Int_t num = atoi(ptr);
 246     if (num>=maxHists) {
 247       Error("LoadHistos()", "Found object with too high index: %s %d",name.Data(),maxHists);
 248       gSystem->Exit(123);
 249     }
 250     //obj->SetDirectory(0);
 251     if (name[1]=='S')
 252       sHists[num]=obj;
 253     else if (name[1]=='M')
 254       mHists[num]=obj;
 255     else {
 256       Error("LoadHistos()", "Found object with unknown: %s %d",name.Data(),maxHists);
 257       gSystem->Exit(123);
 258     }
 259   }
 260
 261   for (int i=0;i<maxHists;++i) {
 262     if (mHists[i]==0)
 263       continue;
 264     if (sHists[i]==0) {
 265       Error("LoadHistos()", "Both histogram ptr should be set: %d",i);
 266       continue;
 267     }
 268
 269     int nEmptyBinsS = 0;
 270     int nEmptyBinsM = 0;
 271     for (int nx=1; nx<36; nx++) {
 272       for (int ny=1; ny<20; ny++) {
 273         if(sHists[i]->GetBinContent(nx,ny)==0)
 274           nEmptyBinsS++;
 275         if(mHists[i]->GetBinContent(nx,ny)==0)
 276           nEmptyBinsM++;
 277       }
 278     }
 279     double emptyFracS = double(nEmptyBinsS)/720.;
 280     double emptyFracM = double(nEmptyBinsM)/720.;
 281     if (emptyFracS > 0.50)
 282       Warning("LoadHistos","%d %s: %.0f%% of bins are EMPTY", i, sHists[i]->GetTitle(), 100*emptyFracS);
 283     if (emptyFracM > 0.50)
 284       Warning("LoadHistos","%d %s: %.0f%% of bins are EMPTY", i, mHists[i]->GetTitle(), 100*emptyFracM);
 285
 286     cHists[i] = (TH2F*)sHists[i]->Clone(Form("hC%d",i));
 287     cHists[i]->Divide(sHists[i], mHists[i], 1./sHists[i]->Integral(), 1./mHists[i]->Integral());
 288   }
 289 }
 290
 291 TH2 *GetHist(int i, int j, int z, int k, char type)
 292 {
 293   if (!fIndex)
 294     return 0;
 295
 296   Int_t ind = fIndex->Eval(i,j,z,k);
 297   TH2* h = 0;
 298
 299   if ((ind>=maxHists)||(ind<0)) {
 300     Error("GetHist", "Input %d %d %d %d gives to large index %d",i,j,z,k,ind);
 301     return 0;
 302   }
 303
 304   if (type == 's')
 305     h = sHists[ind];
 306   else if (type == 'm')
 307     h = mHists[ind];
 308   else if (type == 'c')
 309     h = cHists[ind];
 310   else {
 311     Error("GetHist", "Unknown type %c",type);
 312     return 0;
 313   }
 314
 315   if (h)
 316     if (h->GetSumw2N() == 0)
 317       h->Sumw2();
 318
 319   return h;
 320 }
 321
 322 void SumAndProjectTH2s(Double_t dEtaInnerArg = 1.2, Double_t dEtaOuterArg = 1.799)
 323 {
 324   if (ispp) {
 325     nCentBins = 1;
 326     cb1[0] = 0;
 327     cb2[0] = 0;
 328   }
 329
 330   dEtaInner = dEtaInnerArg;
 331   dEtaOuter = dEtaOuterArg;
 332
 333   outFile = new TFile(Form("$MYDATA/lhc11a_etamin%02d.root", (int)(10*dEtaInner)), "recreate");
 334
 335   LoadHistos();
 336   MakeWeightHistos();
 337   PrintBins();
 338
 339   int nk = ispp ? 1 : 10;
 340   for (int k=0; k<nk; k++) {
 341     TH2* h = GetHist(0,0,4, k, 's');
 342     cout << Form("%s %.2g",h->GetTitle(), h->Integral()) << endl;
 343   }
 344
 345   // Dude, this is such a hack. TODO: calculate safely the number of y
 346   // bins included in an x projection, or vice versa.
 347   double nProjectedBins[] = {8, 12, 18, 18, 18};
 348
 349   kPhiMin[NSJET] = -TMath::PiOver2();
 350   kPhiMax[NSJET] = 3*TMath::PiOver2();
 351   kEtaMin[NSJET] = -dEtaInner;
 352   kEtaMax[NSJET] = +dEtaInner;
 353   kPhiMin[RIDGE] = -TMath::PiOver2();
 354   kPhiMax[RIDGE] = 3*TMath::PiOver2();
 355   kEtaMin[RIDGE] = -dEtaInner; // selection gets inverted
 356   kEtaMax[RIDGE] = +dEtaInner;
 357   kPhiMin[ALL]   = -TMath::PiOver2();
 358   kPhiMax[ALL]   = 3*TMath::PiOver2();
 359   kEtaMin[ALL]   = -dEtaOuter;
 360   kEtaMax[ALL]   = +dEtaOuter;
 361
 362   kPhiMin[ETA_NS]   = -TMath::PiOver2();
 363   kPhiMax[ETA_NS]   = +TMath::PiOver2();
 364   kEtaMin[ETA_NS]   = -dEtaOuter;
 365   kEtaMax[ETA_NS]   = +dEtaOuter;
 366
 367   kPhiMin[ETA_AS]   =   TMath::PiOver2();
 368   kPhiMax[ETA_AS]   = 3*TMath::PiOver2();
 369   kEtaMin[ETA_AS]   = -dEtaOuter;
 370   kEtaMax[ETA_AS]   = +dEtaOuter;
 371
 372
 373   // First loop: sum z-vertex bins - 2D
 374   for (int cb=0; cb<nCentBinsHist; cb++) {
 375
 376     double z_int = hz[cb]->Integral();
 377     double centLo = cbh1[cb];
 378     double centHi = cbh2[cb];
 379     if (ispp) {
 380       centLo = -1;
 381       centHi = 101;
 382     }
 383     cout << Form("Cent bin %d/%d: %.2g to %.2g%%", cb, nCentBinsHist, centLo, centHi)
 384          << endl;
 385
 386     for (int i=0; i<gi->GetN(); i++) {
 387       for (int j=0; j<=i; j++) {
 388
 389         for (int ict=0; ict<nCorrTypes-1; ict++) {
 390
 391           TH2* hz2 = 0;
 392           double wtSum = 0;
 393           bool printZWeights = (0 && i==0 && j==0 && ict==kCA);
 394
 395           if (printZWeights)
 396             cout << Form("i%d j%d %.0fto%.0f: ", i, j, centLo, centHi);
 397
 398           for (int iz=0; iz<nZbins; iz++) {
 399
 400             int zbin1 = hz[cb]->FindBin(z1[iz]);
 401             int zbin2 = hz[cb]->FindBin(z2[iz]-0.001);
 402
 403             double weight = hz[cb]->Integral(zbin1, zbin2) / z_int;
 404
 405             if (ict==kS || ict==kM)
 406               weight = 1.0; // Only correlation functions should be weighted
 407
 408             wtSum += weight;
 409
 410             // There was an off-by-1 bug here!! Getting iz offset from
 411             // zero, but loop was from iz=1 up. Fixed.
 412             TH2* hc2 = GetHist(i,j,iz,cb,sCorrType[ict][0]);
 413             if (!hc2) {
 414               if (iz==0 && i==0 && j==0)
 415                 Info("Loop 1", "Histo not found - skipping");
 416               continue;
 417             }
 418             // if (i==7 && j==6 && ict < kCA) {
 419             //   cout << Form("%s ===>  %.3g", hc2->GetTitle(), hc2->Integral()) << endl;
 420             // }
 421             if (printZWeights) {
 422               cout << Form("%.2gto%.2g:%.2g ",z1[iz], z2[iz], weight);
 423             }
 424
 425             int cLo = ispp ? 0 : centLo;
 426             int cHi = ispp ? 0 : centHi;
 427             const char* newName = Form("%s_%s_%d_%d_%dto%d",
 428                                        "ETAPHI", sCorrType[ict].Data(),
 429                                        i, j, cLo, cHi);
 430
 431             if (iz==0) {
 432               hz2 = (TH2*)hc2->Clone(newName);
 433               hz2->Reset();
 434             }
 435             hz2->Add(hc2, weight);
 436           } // iz loop
 437
 438           if (printZWeights)
 439             cout << endl;
 440
 441           hists->Add(hz2);
 442           if((ict==kCA || ict==kCB) && TMath::Abs(wtSum - 1.) > 0.01)
 443             Warning("MakeWeightHistos()", "weight sum %.2g", wtSum);
 444           if (ict == kCA) {
 445             double normInt = hz2->Integral()/720;
 446             if (TMath::Abs(1-normInt) > 0.05)
 447               Warning("Loop 1", "%s (%d) integral not 1: %.3g",
 448                       hz2->GetName(), (int)fIndex->Eval(i,j,3,cb), normInt);
 449
 450           }
 451
 452         }
 453       }
 454     }
 455   }
 456
 457   cout<<"\nFinished z-vertex bin sums\n"<<endl;
 458
 459   /*
 460   for (int n=0; n<hists->GetEntries(); n++) {
 461     TH2* h = (TH2*)hists->At(n);
 462     TString name(h->GetName());
 463     if (name.Contains("s_0_0"))
 464       cout<<Form("%s %s %.2g",name.Data(), h->GetTitle(), h->Integral())<<endl;
 465   }
 466   */
 467
 468   outFile->mkdir("PbPb");
 469   outFile->cd("PbPb");
 470
 471   // Second loop: sum centrality bins - 2D histos. Purpose: make the
 472   // desired centrality bin k in the TGraph from smaller centrality
 473   // bins cb that already exist.
 474   for (int i=0; i<gi->GetN(); i++) {
 475     for (int j=0; j<=i; j++) {
 476       for (int k=0; k<gk->GetN(); k++) { // the target centrality bin
 477         for (int ict=0; ict<nCorrTypes-1; ict++) {
 478
 479
 480
 481           TObjArray* arr = 0;
 482           std::vector<double> weights;
 483           std::vector<TString> centIntervals;
 484           double wtsum = 0;
 485           // double centLo = cb1[k];
 486           // double centHi = cb2[k]; same as below
 487           double centLo = ispp? 0 : gk->GetX()[k] - gk->GetEX()[k];
 488           double centHi = ispp? 0 : gk->GetX()[k] + gk->GetEX()[k];
 489           const char* name = Form("%s_%s_%d_%d_%.0fto%.0f",
 490                                   "ETAPHI", sCorrType[ict].Data(), i, j, centLo,centHi);
 491           const char* newName = Form("%s_%s_%d_%d_%d", "ETAPHI", sCorrType[ict].Data(), i, j, k);
 492           const char* title   = Form("p_{T}^{t} %s, p_{T}^{a} %s, %.0f-%.0f%%",
 493                                      BinLabel(0, i),
 494                                      BinLabel(1, j),
 495                                      centLo, centHi);
 496
 497           //      Printf("%g %g", centLo, centHi);
 498
 499           TH2* h2 = 0;
 500           h2 = (TH2*) hists->FindObject(name);
 501
 502           if (h2) {
 503             // Give correct overall normalization for CFs.
 504             if (ict == kCA || ict == kCB)
 505               h2->Scale(h2->GetNbinsX()*h2->GetNbinsY()/h2->Integral());
 506
 507             h2->SetNameTitle(newName, title);
 508             h2->Write(newName);
 509
 510             if (1) cout<<h2->GetName()<<" "<<endl;
 511           }
 512           else { // If cent. bin k doesn't exist, it needs to be created as a sum.
 513             for (int cb=0; cb<nCentBinsHist; cb++) { // pre-existing cent bins
 514               const char* name2 = Form("%s_%s_%d_%d_%.0fto%.0f",
 515                                        "ETAPHI", sCorrType[ict].Data(), i, j, cbh1[cb], cbh2[cb]);
 516               TH2* hcb2 = (TH2*) hists->FindObject(name2);
 517               //              TH2* h = GetHist(0,0,4, k, 's');
 518               if (!hcb2) {
 519                 Error("Loop 2","!hcb2");
 520                 continue;
 521               }
 522               cout << name2 <<" "<< cbh1[cb] << " " << centLo << " " << cbh2[cb] << " " << centHi << endl;
 523
 524               // Store relevant histos and their weights
 525               if ( cbh1[cb] >= centLo && cbh2[cb] <= centHi ) {
 526                 if (1) cout<<"Adding " << hcb2->GetName()<<" "<<endl;
 527                 if (!arr) arr = new TObjArray();
 528                 arr->Add(hcb2);
 529
 530                 // weight by integral of same-event 2D distribution rather than event count
 531                 // small but visible correction from old way.
 532                 /* old way:
 533                 weights.push_back(hcent->GetBinContent(cb+1));
 534                 wtsum += hcent->GetBinContent(cb+1);
 535                 */
 536                 // new way
 537                 const char* ss = Form("ETAPHI_s_%d_%d_%.0fto%.0f", i, j, cbh1[cb], cbh2[cb]);
 538                 TH2* hs = (TH2*) hists->FindObject(ss);
 539                 if (!hs)
 540                   Error("Loop 2","Problem finding %s", ss);
 541                 double wt = hs->Integral();
 542                 weights.push_back(wt);
 543                 centIntervals.push_back(TString(Form("%.0fto%.0f", cbh1[cb], cbh2[cb])));
 544                 wtsum += wt;
 545
 546               }
 547             } // cb loop
 548
 549             // Add up histos in arr. This part of code only gets used if
 550             // centrality combination is requested by including new
 551             // bins in cb1,2[] arrays.
 552
 553             bool printWeights = (ict == kCA); // print once/bin, not 3x
 554             double tot=0;
 555
 556             if (printWeights) cout<<name<<" "<<flush;
 557
 558           cout<< arr <<endl;
 559           if (!arr)
 560             continue;
 561
 562             for (int n=0; n<arr->GetEntries(); n++) {
 563
 564               tot += weights.at(n)/wtsum;
 565               TH2* hn = (TH2*)arr->At(n);
 566
 567               if (n==0) {
 568                 h2 = (TH2*)hn->Clone(newName);
 569                 h2->Scale(weights.at(n)/wtsum);
 570               }
 571               else {
 572                 h2->Add(hn, weights.at(n)/wtsum);
 573               }
 574               if (printWeights)
 575                 cout << centIntervals.at(n).Data() << ":"
 576                      << weights.at(n)/wtsum << " " << flush;
 577             }
 578             if (printWeights) cout << " = " << tot << endl;
 579
 580             // Give correct overall normalization for CFs.
 581             if (ict == kCA || ict == kCB)
 582               h2->Scale(h2->GetNbinsX()*h2->GetNbinsY()/h2->Integral());
 583
 584             h2->SetNameTitle(newName, title);
 585             h2->Write(newName);
 586
 587             if (arr)
 588               delete arr; // arr doesn't own its histos, so ok to delete.
 589           }
 590         } // ict
 591       } // k
 592     } // j
 593   } // i
 594
 595   cout<<"\nFinished centrality bin sums\n"<<endl;
 596
 597   outFile->cd();
 598   for (int ijkz=0; ijkz<4; ijkz++)
 599     BinGraph(ijkz)->Write();
 600   profs->Write();
 601
 602   TObjArray* summed = GetHistList(*outFile, "TH2", "PbPb");
 603   TObjArray* dphis = new TObjArray();
 604
 605   outFile->cd("PbPb");
 606   for (int n=0; n<summed->GetEntries(); n++) {
 607     TH2* h = (TH2*)summed->At(n);
 608     for (int r=0; r<5; r++) {
 609       TString name1D(h->GetName());
 610       name1D.ReplaceAll("ETAPHI", Region(r));
 611       TString projStr = "x";
 612       if (r==RIDGE)
 613         projStr.Prepend("-");
 614       else if (r==ETA_NS || r==ETA_AS)
 615         projStr = "y";
 616
 617       TH1* hp = ProjectTH2(*h, projStr,
 618                            kPhiMin[r], kEtaMin[r], kPhiMax[r], kEtaMax[r],
 619                            name1D, "e");
 620       hp->SetTitle(h->GetTitle());
 621
 622       // TODO check this
 623       if (name1D.Contains("_c_"))
 624         hp->Scale(1./nProjectedBins[r]);
 625
 626       if (r==NSJET || r==RIDGE || r==ALL)
 627         dphis->Add(hp);
 628
 629       hp->Write();
 630     }
 631   }
 632
 633   // CF defs A and B
 634   // First find existing "Def. B" CF _c_. Then create "Def. A" CFs as well.
 635   for (int n=0; n<dphis->GetEntries(); n++) {
 636     TObject* obj = dphis->At(n);
 637     TString hname =  obj->GetName();
 638     TH1* hcB = (TH1*)obj;
 639     if (hname.Contains("_c_") ) {
 640       hcB = (TH1*)obj;
 641     }
 642     else
 643       continue;
 644
 645     TString s(hname);
 646     TString m(hname);
 647     TString cA(hname);
 648     TString cB(hname);
 649     s.ReplaceAll("_c_", "_s_");
 650     m.ReplaceAll("_c_", "_m_");
 651     cA.ReplaceAll("_c_", "_cA_");
 652     cB.ReplaceAll("_c_", "_cB_");
 653
 654     TH1* hsame  = (TH1*)dphis->FindObject(s.Data());
 655     TH1* hmixed = (TH1*)dphis->FindObject(m.Data());
 656     TH1* hcA = (TH1*)hcB->Clone(cA.Data());
 657
 658     if (0)
 659       Info("", "%s %s %s", hsame->GetName(), hmixed->GetName(), hcA->GetName());
 660
 661     hcA->Divide(hsame, hmixed, 1./hsame->Integral(), 1./hmixed->Integral());
 662     hcB->SetName(cB.Data());
 663
 664     hcA->Write();
 665     hcB->Write();
 666
 667   }
 668
 669   outFile->cd("PbPb");
 670   int nObjectsInFile = gDirectory->GetListOfKeys()->GetEntries();
 671   cout <<"\nWrote " << nObjectsInFile
 672        << " objects to PbPb/ in " << outFile->GetName() << endl;
 673
 674   cout << "Closing file..." << endl;
 675   outFile->Close();
 676
 677   return;
 678 }
 679
 680 const char* Region(int r)
 681 {
 682   return regionStr[r];
 683 }
 684
 685 char* BinLabel(int ijkz, int i)
 686 {
 687   // Returns bin extents like "a-b" to <= 4 sig. fig. precision
 688   TGraphErrors* g = BinGraph(ijkz);
 689   if (i<0 || i>=g->GetN()) {
 690     Error("bin_label((TGraphErrors&, Int_t)", "Error: no bin %d", i);
 691     return 0;
 692   }
 693
 694   double x1 = g->GetX()[i] - g->GetEX()[i];
 695   double x2 = g->GetX()[i] + g->GetEX()[i];
 696
 697   return Form("%.*g-%.*g", 4, x1, 4, x2);
 698 }
 699
 700 TGraphErrors* BinGraph(int ijkz)
 701 {
 702   if (ijkz < 0 || ijkz > 3)
 703     return 0;
 704   else if (ijkz==0) return gi;
 705   else if (ijkz==1) return gj;
 706   else if (ijkz==2) return gk;
 707   else if (ijkz==3) return gz;
 708   else
 709     return 0;
 710 }
 711
 712 TObjArray* GetHistList(TFile& file, TString clname, TString dir)
 713 {
 714   file.cd(dir.Data());
 715
 716   TObjArray* hList = new TObjArray();
 717   TIter next(gDirectory->GetListOfKeys());
 718   TKey *key;
 719
 720   while ((key=(TKey*)next())) {
 721     TString className(key->GetClassName());
 722     TString keyName(key->GetName());
 723     if (0)
 724       printf("%10s %20s\n", className.Data(), keyName.Data());
 725
 726     if (className.Contains(clname) && clname.Contains("TH2")) {
 727       hList->Add((TH2*)gDirectory->Get(keyName.Data()));
 728     }
 729   }
 730
 731   cout << hList->GetEntries() << " objects retrieved from "
 732        << file.GetName()  << "/" << gDirectory->GetName()
 733        << endl;
 734
 735   return hList;
 736 }
 737
 738 TH1*
 739 ProjectTH2(TH2& h, TString xy, float x1, float y1, float x2, float y2,
 740            TString name, TString opt)
 741 {
 742   TH1* hp = 0;
 743   double x[] = {x1,x2,x2,x1,x1};
 744   double y[] = {y1,y1,y2,y2,y1};
 745   TCutG *cutg = new TCutG("cutg", 5, x, y);
 746   double area = cutg->Area();
 747   TString projOpt = "[cutg]";
 748
 749   // Invert selection: fill bins outside rectangle
 750   if (xy.Contains("-"))
 751     projOpt.ReplaceAll("[cutg]", "[-cutg]");
 752
 753   // "e" errors, "d" draw in current pad, "o" original full axis
 754   projOpt.Append(opt);
 755
 756   if (xy.Contains("x")) {
 757     hp = h.ProjectionX(name.Data(), 0, -1, projOpt.Data());
 758   }
 759   else if (xy.Contains("y")) {
 760     hp = h.ProjectionY(name.Data(), 0, -1, projOpt.Data());
 761   }
 762   if (!hp)
 763     Error("ProjectTH2()","!hp");
 764
 765
 766   if (0) { // For test/debug purposes
 767     int bx1, bx2, by1, by2;
 768     double wx, hy;
 769     TAxis *ax = h.GetXaxis(), *ay = h.GetYaxis();
 770
 771     bx1 = ax->FindBin(x1);
 772     bx2 = ax->FindBin(x2);
 773     by1 = ay->FindBin(y1);
 774     by2 = ay->FindBin(y2);
 775     wx = ax->GetBinUpEdge(bx2) - ax->GetBinLowEdge(bx1);
 776     hy = ay->GetBinUpEdge(by2) - ay->GetBinLowEdge(by1);
 777
 778     if (xy.Contains("-"))
 779       cout << Form("Selection is inverted.") << endl;
 780
 781     cout << Form("Selection area: %f", area) << endl;
 782     cout << Form("Actual projected area %f", wx*hy) << endl;
 783   }
 784
 785   return hp;
 786 }
 787
 788 void PrintBins()
 789 {
 790   TGraphErrors* g = 0;
 791
 792   cout << endl;
 793   cout <<"-----------------Bin Information-----------------"<< endl;
 794   for (int n=0; n<4; n++) {
 795     g = BinGraph(n);
 796     cout << g->GetName() << ": ";
 797     for (int i=0; i<g->GetN(); i++) {
 798       cout << Form("%s ", BinLabel(n, i));
 799     }
 800     cout << endl;
 801   }
 802   cout <<"-------------------------------------------------"<< endl;
 803  return;
 804 }