AliPhysics  6bc8652 (6bc8652)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SimpleCorrect.C
Go to the documentation of this file.
1 #ifndef __CINT__
2 # include <TCanvas.h>
3 # include <TList.h>
4 # include <TH1.h>
5 # include <TH2.h>
6 # include <THStack.h>
7 # include <TFile.h>
8 # include <TError.h>
9 # include <TLatex.h>
10 # include <TFitResult.h>
11 # include <TF1.h>
12 # include <TMath.h>
13 # include <THashList.h>
14 #else
15 class TList;
16 class TSeqCollection;
17 class TH1;
18 class TH2;
19 class THStack;
20 class TAxis;
21 class TDirectory;
22 #endif
23 
24 class TCanvas; // Auto load
25 
27 
31 const char* dndetaTitle = "d#it{N}_{ch}/d#it{#eta}";
35 const Color_t cc[] = { kMagenta+2,
36  kBlue+2,
37  kAzure-1, // 10,
38  kCyan+2,
39  kGreen+1,
40  kSpring+5,//+10,
41  kYellow+1,
42  kOrange+5,//+10,
43  kRed+1,
44  kPink+5,//+10,
45  kBlack };
46 
47 //====================================================================
61 TObject* GetO(TSeqCollection* l, const char* name, TClass* cls)
62 {
63  if (!l) {
64  Warning("GetO", "No list passed");
65  return 0;
66  }
67  TObject* o = l->FindObject(name);
68  if (!o) {
69  Warning("GetO", "No object named %s found in list %s",
70  name, l->GetName());
71  return 0;
72  }
73  if (!cls) return o;
74 
75  if (!o->IsA()->InheritsFrom(cls)) {
76  Warning("GetO", "Object %s read from %s is not a %s but a %s",
77  name, l->GetName(), cls->GetName(), o->ClassName());
78  return 0;
79  }
80  return o;
81 }
82 //____________________________________________________________________
83 TObject* CopyO(TSeqCollection* l, const char* name, const char* newName,
84  TClass* cls)
85 {
86  TObject* orig = GetO(l, name, cls);
87  if (!orig) return 0;
88  TObject* copy = orig->Clone(newName ? newName : name);
89  return copy;
90 }
91 //____________________________________________________________________
100 TH1* GetH1(TSeqCollection* l, const char* name)
101 {
102  return static_cast<TH1*>(GetO(l,name,TH1::Class()));
103 }
104 //____________________________________________________________________
113 TH2* GetH2(TSeqCollection* l, const char* name)
114 {
115  return static_cast<TH2*>(GetO(l,name,TH2::Class()));
116 }
117 //____________________________________________________________________
127 TH1* CopyH1(TSeqCollection* l, const char* name, const char* newName=0)
128 {
129  TH1* copy = static_cast<TH1*>(CopyO(l,name,newName,TH1::Class()));
130  if (copy) copy->SetDirectory(0);
131  return copy;
132 }
133 //____________________________________________________________________
143 TH2* CopyH2(TSeqCollection* l, const char* name, const char* newName=0)
144 {
145  TH2* copy = static_cast<TH2*>(CopyO(l,name,newName,TH2::Class()));
146  if (copy) copy->SetDirectory(0);
147  return copy;
148 }
149 /* @} */
150 
151 //====================================================================
165 Bool_t CheckAxisNBins(const char* which,
166  const TAxis* a1,
167  const TAxis* a2)
168 {
169  if (a1->GetNbins() != a2->GetNbins()) {
170  ::Warning("CheckAxisNBins", "Incompatible number %s bins: %d vs %d",
171  which, a1->GetNbins(), a2->GetNbins());
172  return false;
173  }
174  return true;
175 }
176 //____________________________________________________________________
186 Bool_t CheckAxisLimits(const char* which,
187  const TAxis* a1,
188  const TAxis* a2)
189 {
190  if (!TMath::AreEqualRel(a1->GetXmin(), a2->GetXmin(),1.E-12) ||
191  !TMath::AreEqualRel(a1->GetXmax(), a2->GetXmax(),1.E-12)) {
192  Warning("CheckAxisLimits",
193  "Limits of %s axis incompatible [%f,%f] vs [%f,%f]", which,
194  a1->GetXmin(), a1->GetXmax(), a2->GetXmin(), a2->GetXmax());
195  return false;
196  }
197  return true;
198 }
199 //____________________________________________________________________
209 Bool_t CheckAxisBins(const char* which,
210  const TAxis* a1,
211  const TAxis* a2)
212 {
213  const TArrayD * h1Array = a1->GetXbins();
214  const TArrayD * h2Array = a2->GetXbins();
215  Int_t fN = h1Array->fN;
216  if ( fN == 0 ) return true;
217  if (h2Array->fN != fN) {
218  // Redundant?
219  Warning("CheckAxisBins", "Not equal number of %s bin limits: %d vs %d",
220  which, fN, h2Array->fN);
221  return false;
222  }
223  else {
224  for (int i = 0; i < fN; ++i) {
225  if (!TMath::AreEqualRel(h1Array->GetAt(i),h2Array->GetAt(i),1E-10)) {
226  Warning("CheckAxisBins",
227  "%s limit # %3d incompatible: %f vs %f",
228  which, i, h1Array->GetAt(i),h2Array->GetAt(i));
229  return false;
230  }
231  }
232  }
233  return true;
234 }
235 //____________________________________________________________________
247 Bool_t CheckAxisLabels(const char* which,
248  const TAxis* a1,
249  const TAxis* a2)
250 {
251  // check that axis have same labels
252  THashList *l1 = (const_cast<TAxis*>(a1))->GetLabels();
253  THashList *l2 = (const_cast<TAxis*>(a2))->GetLabels();
254 
255  if (!l1 && !l2) return true;
256  if (!l1 || !l2) {
257  Warning("CheckAxisLabels", "Difference in %s labels: %p vs %p",
258  which, l1, l2);
259  return false;
260  }
261  // check now labels sizes are the same
262  if (l1->GetSize() != l2->GetSize()) {
263  Warning("CheckAxisLabels", "Different number of %s labels: %d vs %d",
264  which, l1->GetSize(), l2->GetSize());
265  return false;
266  }
267  for (int i = 1; i <= a1->GetNbins(); ++i) {
268  TString label1 = a1->GetBinLabel(i);
269  TString label2 = a2->GetBinLabel(i);
270  if (label1 != label2) {
271  Warning("CheckAxisLabels", "%s label # %d not the same: '%s' vs '%s'",
272  which, i, label1.Data(), label2.Data());
273  return false;
274  }
275  }
276 
277  return true;
278 }
279 //____________________________________________________________________
291 Bool_t CheckAxis(const char* which,
292  const TAxis* a1,
293  const TAxis* a2,
294  Bool_t alsoLbls)
295 {
296  if (!CheckAxisNBins (which, a1, a2)) return false;
297  if (!CheckAxisLimits(which, a1, a2)) return false;
298  if (!CheckAxisBins (which, a1, a2)) return false;
299  if (alsoLbls && !CheckAxisLabels(which, a1, a2)) return false;
300  return true;
301 }
302 //____________________________________________________________________
311 Bool_t CheckConsistency(const TH1* h1, const TH1* h2)
312 {
313  // Check histogram compatibility
314  if (h1 == h2) return true;
315 
316  if (h1->GetDimension() != h2->GetDimension()) {
317  Warning("CheckConsistency",
318  "%s and %s have different dimensions %d vs %d",
319  h1->GetName(), h2->GetName(),
320  h1->GetDimension(), h2->GetDimension());
321  return false;
322  }
323  Int_t dim = h1->GetDimension();
324 
325  Bool_t ret = true;
326  Bool_t alsoLbls = (h1->GetEntries() != 0 && h2->GetEntries() != 0);
327  if (!CheckAxis("X", h1->GetXaxis(), h2->GetXaxis(), alsoLbls)) ret = false;
328  if (dim > 1 &&
329  !CheckAxis("Y", h1->GetYaxis(), h2->GetYaxis(), alsoLbls)) ret = false;
330  if (dim > 2 &&
331  !CheckAxis("Z", h1->GetZaxis(), h2->GetZaxis(), alsoLbls)) ret = false;
332 
333  return ret;
334 }
335 
336 //____________________________________________________________________
345 TH1* Coerce(const TH1* templ, TH1* src)
346 {
347  if (CheckConsistency(templ, src)) return src;
348  TH1* ret = static_cast<TH1*>(templ->Clone(src->GetName()));
349  ret->SetDirectory(0);
350  ret->Reset();
351  for (Int_t i = 1; i <= ret->GetNbinsX(); i++) {
352  Double_t x = ret->GetXaxis()->GetBinCenter(i);
353  Int_t b = src->GetXaxis()->FindBin(x);
354  if (b <= 0 || b >= src->GetNbinsX()) continue;
355  Double_t c = src->GetBinContent(b);
356  Double_t e = src->GetBinError(b);
357  ret->SetBinContent(i, c);
358  ret->SetBinError (i, e);
359  }
360  return ret;
361 }
362 
363 //____________________________________________________________________
372 TH2* Coerce(const TH2* templ, TH2* src)
373 {
374  if (CheckConsistency(templ, src)) return src;
375  TH2* ret = static_cast<TH2*>(templ->Clone(src->GetName()));
376  ret->SetDirectory(0);
377  ret->Reset();
378  for (Int_t i = 1; i <= ret->GetNbinsX(); i++) {
379  Double_t x = ret->GetXaxis()->GetBinCenter(i);
380  Int_t bx = src->GetXaxis()->FindBin(x);
381  if (bx <= 0 || bx >= src->GetNbinsX()) continue;
382  for (Int_t j = 1; j <= ret->GetNbinsY(); j++) {
383  Double_t y = ret->GetYaxis()->GetBinCenter(j);
384  Int_t by = src->GetYaxis()->FindBin(y);
385  Double_t c = src->GetBinContent(bx,by);
386  Double_t e = src->GetBinError(bx,by);
387  ret->SetBinContent(i, j, c);
388  ret->SetBinError (i, j, e);
389  }
390  }
391  return ret;
392 }
393 /* @} */
394 
395 //====================================================================
414  Color_t color,
415  Style_t marker=20,
416  Double_t size=1.,
417  Style_t fill=0,
418  Style_t line=1,
419  Width_t width=1)
420 {
421  if (!h) return 0;
422  h->SetMarkerColor(color);
423  h->SetMarkerStyle(marker);
424  h->SetMarkerSize(size);
425  h->SetFillColor(color);
426  h->SetFillStyle(fill);
427  h->SetLineColor(color);
428  h->SetLineStyle(line);
429  h->SetLineWidth(width);
430  h->GetXaxis()->SetNdivisions(210);
431  h->GetYaxis()->SetNdivisions(210);
432 }
433 
434 //____________________________________________________________________
441 void PrintH(TH2* h, Int_t prec=2)
442 {
443  Printf("Content of %s - %s", h->GetName(), h->GetTitle());
444  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
445  printf("%3d: ", i);
446  for (Int_t j = 1; j <= h->GetNbinsY(); j++) {
447  Double_t c = h->GetBinContent(i,j);
448  Double_t e = h->GetBinError (i,j);
449  if (TMath::IsNaN(c) || TMath::IsNaN(e))
450  printf("*** NAN ***");
451  else
452  printf("%.*f+/-%.*f ", prec, c, prec, e);
453  }
454  printf("\n");
455  }
456 }
457 //____________________________________________________________________
464 void PrintH(TH1* h, Int_t prec=2)
465 {
466  Printf("Content of %s - %s", h->GetName(), h->GetTitle());
467  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
468  if (h->GetBinContent(i) <= 1e-6) continue;
469  printf("%3d (%+6.3f): %.*f+/-%.*f\n", i,
470  h->GetXaxis()->GetBinCenter(i),
471  prec, h->GetBinContent(i),
472  prec, h->GetBinError(i));
473  }
474 }
475 /* @} */
476 //====================================================================
496  Double_t min,
497  Double_t max,
498  Double_t& err)
499 {
500  const Double_t epsilon = 1e-6;
501  Int_t bMin = h->FindBin(min+epsilon);
502  Int_t bMax = h->FindBin(max-epsilon);
503  if (bMin < 1) bMin = 0; // Include underflow bin
504  Double_t val = h->IntegralAndError(bMin, bMax, err);
505  // For a-symmetric histograms, return
506  if (TMath::Abs(h->GetXaxis()->GetXmin()+h->GetXaxis()->GetXmax())>=epsilon)
507  return val;
508 
509  // Otherwise, also integrate negative side
510  Double_t err2;
511  bMin = h->FindBin(-min+epsilon);
512  bMax = h->FindBin(-max-epsilon);
513  val += h->IntegralAndError(bMin, bMax, err2);
514  err = TMath::Sqrt(err*err+err2*err2);
515  return val;
516 }
517 //____________________________________________________________________
534  Double_t d, Double_t ed,
535  Double_t& er)
536 {
537  Double_t r = 0;
538  er = 0;
539  if (TMath::Abs(n) < 1e-16 || TMath::Abs(d) < 1e-16) return 0;
540  r = n/d;
541  er = TMath::Sqrt(en*en/n/n + ed*ed/d/d);
542  return r;
543 }
544 /* @} */
545 
546 //====================================================================
551 //____________________________________________________________________
567 TH2* Scale(TH2* h, TH1* g)
568 {
569  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
570  Double_t dc = g->GetBinContent(i);
571  Double_t de = g->GetBinError (i);
572  Double_t dr = (dc > 1e-6 ? 1/dc : 0);
573  Double_t dq = dr * de;
574  for (Int_t j = 1; j <= h->GetNbinsY(); j++) {
575  Double_t nc = h->GetBinContent(i,j);
576  Double_t ne = h->GetBinError (i,j);
577  Double_t ns = (nc > 0 ? ne/nc : 0);
578  Double_t sc = dc * nc;
579  Double_t se = sc*TMath::Sqrt(ns*ns+dq*dq);
580  // Printf("Setting bin %2d,%2d to %f +/- %f", sc, se);
581  h->SetBinContent(i,j,sc);
582  h->SetBinError (i,j,se);
583  }
584  }
585  return h;
586 }
587 //____________________________________________________________________
605 {
606  Double_t rr = se/s;
607  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
608  Double_t c = h->GetBinContent(i);
609  Double_t e = h->GetBinError (i);
610  Double_t es = (c > 0 ? e/c : 0);
611  Double_t sc = s * c;
612  Double_t se = sc*TMath::Sqrt(es*es+rr*rr);
613  // Printf("Setting bin %2d,%2d to %f +/- %f", sc, se);
614  h->SetBinContent(i,sc);
615  h->SetBinError (i,se);
616  }
617  return h;
618 }
619 
620 //____________________________________________________________________
638 {
639  Double_t rr = se/s;
640  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
641  for (Int_t j = 1; j <= h->GetNbinsY(); j++) {
642  Double_t c = h->GetBinContent(i,j);
643  Double_t e = h->GetBinError (i,j);
644  Double_t es = (c > 0 ? e/c : 0);
645  Double_t sc = s * c;
646  Double_t se = sc*TMath::Sqrt(es*es+rr*rr);
647  // Printf("Setting bin %2d,%2d to %f +/- %f", sc, se);
648  h->SetBinContent(i,j,sc);
649  h->SetBinError (i,j,se);
650  }
651  }
652  return h;
653 }
664 TH2* CorrectForIPz(TH2* in, TH1* ipz, Bool_t full=true)
665 {
666  // Printf("Scaling to IP of %s", in->GetName());
667  for (Int_t iy = 1; iy <= in->GetNbinsY(); iy++) {
668  Double_t z = in->GetYaxis()->GetBinCenter(iy);
669  Int_t bin = ipz->GetXaxis()->FindBin(z);
670  Double_t nEv = ipz->GetBinContent(bin);
671  Double_t eEv = ipz->GetBinError (bin);
672  Double_t esc = (nEv > 0 ? 1./nEv : 0);
673  Double_t rE2 = esc*esc*eEv*eEv;
674  // Printf("z=%f -> %f +/- %f events", z, nEv, eEv);
675  for (Int_t ix = 1; ix <= in->GetNbinsX(); ix++) {
676  Double_t c = in->GetBinContent(ix,iy);
677  Double_t e = in->GetBinError (ix,iy);
678  Double_t r = (c > 0 ? e/c : 0);
679  // Scale by number of events, and error propagate
680  Double_t sc = c * esc;
681  Double_t se = 0;
682  if (full)se = sc * TMath::Sqrt(r*r+rE2);
683  else se = e * esc;
684  Double_t scl = 1 / in->GetXaxis()->GetBinWidth(ix);
685  // Printf("Setting bin (%3d,%3d) to %f +/- %f (%f %5.1f%% %5.1f%%)",
686  // ix,iy,scl*sc, scl*se, scl, 100*r, 100*esc*eEv);
687  in->SetBinContent(ix, iy, scl*sc);
688  in->SetBinError (ix, iy, scl*se);
689  }
690  }
691  return in;
692 }
693 
694 //____________________________________________________________________
708 TH1* Avg(TH2* h,
709  TH1* ipz,
710  UShort_t mode,
711  const char* name,
712  TH2* other=0)
713 {
714  if (!h) return 0;
715  TH2* mask = other ? other : h;
716  // Printf("Averaging %s - %s (%p) %s", h->GetName(), h->GetTitle(),
717  // ipz, (mode ? "errors only" : "full scale"));
718  Int_t nIPz = h->GetNbinsY();
719  Int_t nEta = h->GetNbinsX();
720  TH1* p = h->ProjectionX(name,1,nIPz,"e");
721  p->SetDirectory(0);
722  p->SetFillColor(0);
723  p->SetFillStyle(0);
724  p->SetYTitle(Form("#LT(%s)#GT", h->GetYaxis()->GetTitle()));
725  p->Reset();
726  for (Int_t etaBin = 1; etaBin <= nEta; etaBin++) {
727  TArrayD hv(nIPz);
728  TArrayD he(nIPz);
729  TArrayD hr(nIPz);
730  TArrayI hb(nIPz);
731  Int_t j = 0;
732  for (Int_t ipzBin = 1; ipzBin <= nIPz; ipzBin++) {
733  Double_t bc = mask->GetBinContent(etaBin, ipzBin);
734  if (bc < 1e-9) continue; // Low value
735  Double_t be = mask->GetBinError (etaBin, ipzBin);
736  if (TMath::IsNaN(be)) continue; // Bad error value
737  hv[j] = bc;
738  he[j] = be;
739  hr[j] = be/bc;
740  hb[j] = ipzBin;
741  j++;
742  }
743  // Now we have all interesting values. Sort the relative error
744  // array to get the most significant first.
745  TArrayI idx(nIPz);
746  TMath::Sort(j, hr.fArray, idx.fArray, false);
747  Double_t nsum = 0; // Running weighted sum
748  Double_t nsume = 0; // Running weighted sum error
749  Double_t dsum = 0;
750  Double_t dsume = 0;
751  Int_t n = 0;
752  Double_t rat = 1e99;
753 
754  Int_t k = 0;
755  for (k = 0; k < j; k++) {
756  Int_t l = idx[k]; // Sorted index
757  Int_t ipzBin = hb[l];
758  Double_t hvv = hv[l];
759  Double_t hee = he[l];
760  Double_t hrr = hr[l];
761  Double_t x = TMath::Sqrt(nsume+hee*hee)/(nsum+hvv);
762  if (x > rat) {
763  continue; // Ignore - does not help
764  }
765  Double_t by = h->GetYaxis()->GetBinCenter(ipzBin);
766  Int_t ib = ipz ? ipz->FindBin(by) : 0;
767  rat = x;
768  nsum += h->GetBinContent(etaBin, ipzBin);
769  nsume += TMath::Power(h->GetBinError(etaBin, ipzBin), 2);
770  // If we do not have the vertex distribution, then just count
771  // number of observations.
772  dsum += !ipz ? 1 : ipz->GetBinContent(ib);
773  dsume += !ipz ? 0 : TMath::Power(ipz->GetBinError(ib),2);
774  n++;
775  }
776  if (k == 0 || n == 0) {
777  // ::Warning("Average", "Eta bin # %3d has no data",etaBin);
778  continue; // This eta bin empty!
779  }
780 
781  Double_t norm = (mode > 0 ? n : dsum);
782  Double_t rne = nsume/nsum/nsum;
783  Double_t rde = dsume/dsum/dsum;
784  Double_t av = nsum/norm;
785  Double_t ave = 0;
786  if (mode==2) ave = TMath::Sqrt(nsume)/n;
787  else ave = av*TMath::Sqrt(rne+rde);
788 #if 0
789  Printf("%10s - bin %3d (%+6.3f) count=%2d n=%9.3f+/-%9.3f d=%9.3f+/-%9.3f "
790  "norm=%9.3f -> %9.3f +/- %9.3f",
791  h->GetName(), etaBin, h->GetXaxis()->GetBinCenter(etaBin), n,
792  nsum, TMath::Sqrt(nsume)/(mode == 2 ? n : nsum),
793  dsum, TMath::Sqrt(dsume)/(mode == 2 ? n : dsum),
794  norm, av, ave);
795 #endif
796  //Info("AverageSim", "Setting eta bin # %3d to %f +/- %f", etaBin,av,ave);
797  p->SetBinContent(etaBin, av);
798  p->SetBinError (etaBin, ave);
799  }
800  if (mode == 0) p->Scale(1, "width");
801  return p;
802 }
803 //____________________________________________________________________
818 TH1* Avg2(TH2* h, TH1* ipz, UShort_t mode, const char* name)
819 {
820  TH2* tmp = static_cast<TH2*>(h->Clone("tmp"));
821  tmp->SetDirectory(0);
822  CorrectForIPz(tmp, ipz, true);
823 
824  TH1* ret = Avg(tmp, 0, 2, name);
825  if (tmp) delete tmp;
826 
827  return ret;
828 }
829 /* @} */
830 
831 //====================================================================
859  TH1* deltaBg,
860  Double_t& eR,
861  Double_t nEvRec=1,
862  Double_t nEvBg =1,
863  Double_t lim=1e9)
864 {
865  if (!deltaRec || !deltaBg) {
866  Error("GetDeltaScale", "Missing input histogram(s): %p %p",
867  deltaRec, deltaBg);
868  return 0;
869  }
870  deltaRec->Scale(1./nEvRec);
871  deltaBg ->Scale(1./nEvBg);
872  Double_t top = TMath::Min(lim,deltaRec->GetXaxis()->GetXmax());
873  Double_t eRec, eBg;
874  Double_t iRec = Integrate(deltaRec, 5, top, eRec);
875  Double_t iBg = Integrate(deltaBg, 5, top, eBg);
876  //Printf("Integral of reconstructed Delta: %14.1f +/- %14.2f", iRec, eRec);
877  //Printf("Integral of injected Delta: %14.1f +/- %14.2f", iBg, eBg);
878  Double_t r = RatioE(iRec, eRec, iBg, eBg, eR);
879 #if 0
880  Printf("%20s: data=%9g+/-%9g inj=%9g+/-%9g -> %9g+/-%9g",
881  list->GetName(), iRec, eRec, iBg, eBg, r, eR);
882 #endif
883  return r;
884 }
885 
886 //____________________________________________________________________
899 Double_t GetDeltaScale(TSeqCollection* list,
900  Int_t b,
901  Double_t& eR,
902  Double_t nEv=1,
903  Double_t lim=1e9)
904 {
905  TH1* deltaRec = CopyH1(list,Form("b%d_TrData_WDist", b));
906  TH1* deltaBg = CopyH1(list,Form("b%d_TrInj_WDist", b));
907  Double_t scale = GetDeltaScale(deltaRec, deltaBg, eR, nEv, nEv,lim);
908  delete deltaRec;
909  delete deltaBg;
910  return scale;
911 }
912 //____________________________________________________________________
928 Double_t GetDeltaScale(TSeqCollection* realList,
929  TSeqCollection* simList,
930  Int_t b,
931  Double_t& eR,
932  Double_t nEvReal=1,
933  Double_t nEvSim=1,
934  Double_t lim=1e9)
935 {
936  TH1* deltaReal = CopyH1(realList,Form("b%d_TrData_WDist", b));
937  TH1* deltaSim = CopyH1(simList, Form("b%d_TrData_WDist", b));
938  Double_t scale = GetDeltaScale(deltaReal,deltaSim,eR,nEvReal,nEvSim,lim);
939  delete deltaReal;
940  delete deltaSim;
941  return scale;
942 }
943 
944 //____________________________________________________________________
958 TH1* GetDeltaDist(TH2* deltaReal,
959  TH2* deltaSim,
960  Double_t nEvReal,
961  Double_t nEvSim,
962  Double_t lim=1e9)
963 {
964  Double_t top = deltaReal->GetYaxis()->GetXmax();
965  TH1* ret = deltaReal->ProjectionX("scalar");
966  ret->Reset();
967  ret->SetTitle("Background scale (#it{k})");
968  ret->SetYTitle("#delta scale");
969  ret->SetXTitle("#it{#eta}");
970  ret->SetDirectory(0);
971  for (Int_t i = 1; i <= ret->GetNbinsX(); i++) {
972  TH1* tmpReal = deltaReal->ProjectionY("tmpReal",i,i);
973  TH1* tmpSim = deltaSim ->ProjectionY("tmpSim", i,i);
974  Double_t eR, r = GetDeltaScale(tmpReal, tmpSim, eR, nEvReal, nEvSim, lim);
975  ret->SetBinContent(i, r);
976  ret->SetBinError (i, eR);
977  delete tmpReal;
978  delete tmpSim;
979  }
980  return ret;
981 }
982 
983 //____________________________________________________________________
995 TH1* GetDeltaDist(TSeqCollection* list,
996  Int_t b,
997  Double_t nEv=1)
998 {
999  TH2* deltaRec = GetH2(list,Form("b%d_TrData_WDvsEta", b));
1000  TH2* deltaBg = GetH2(list,Form("b%d_TrInj_WDvsEta", b));
1001  return GetDeltaDist(deltaRec, deltaBg, nEv, nEv);
1002 }
1003 
1004 //____________________________________________________________________
1019 TH1* GetDeltaDist(TSeqCollection* realList,
1020  TSeqCollection* simList,
1021  Int_t b,
1022  Double_t nEvReal,
1023  Double_t nEvSim,
1024  Double_t lim=1e9)
1025 {
1026  TH2* deltaReal = CopyH2(realList,Form("b%d_TrData_WDvsEta", b));
1027  TH2* deltaSim = CopyH2(simList, Form("b%d_TrData_WDvsEta", b));
1028  return GetDeltaDist(deltaReal, deltaSim, nEvReal, nEvSim, lim);
1029 }
1030 /* @} */
1031 
1032 
1033 
1034 //====================================================================
1053  TSeqCollection* list,
1054  Int_t b,
1055  Double_t nEv=1,
1056  Double_t lim=1e9)
1057 {
1058  Double_t eR = 0;
1059  TH2* in = GetH2(list,Form("b%d_Tr%s_ZvEtaCutT",b,comb ? "Comb" : "Inj"));
1060  TH2* ret = static_cast<TH2*>(in->Clone("bg"));
1061  in ->SetDirectory(0);
1062  ret->SetDirectory(0);
1063  ret->SetStats(0);
1064  if (comb) return ret;
1065 
1066 #if 1
1067  // We should do the below, but that doesn't work as it creates
1068  // negative signals!
1069  Double_t r = GetDeltaScale(list, b, eR, nEv, lim);
1070  // Printf("Ratio of integrals: %14.3f +/- %f14.3", r, eR);
1071  return Scale(ret, r, eR);
1072 #else
1073  TH1* r = GetDeltaDist(list, b, nEv);
1074  Scale(ret, r);
1075  // PrintH(r);
1076  delete r;
1077  // PrintH(ret);
1078  return ret;
1079 #endif
1080 }
1081 //____________________________________________________________________
1103 TH2* CorrectSignal(Bool_t comb, TH2* meas, TH2* bg, TH2* simMeas)
1104 {
1105  if (!comb) {
1106  meas->Add(bg, -1);
1107  return meas;
1108  }
1109  TH2* beta = static_cast<TH2*>(bg->Clone("beta"));
1110  beta->Divide(simMeas);
1111  beta->SetDirectory(0);
1112  TH2* one = static_cast<TH2*>(bg->Clone("one"));
1113  one->Reset();
1114  one->SetDirectory(0);
1115  for (Int_t i = 1; i <= one->GetNbinsX(); i++)
1116  for (Int_t j = 1; j <= one->GetNbinsY(); j++)
1117  one->SetBinContent(i,j,1);
1118  one->Add(beta,-1);
1119  meas->Multiply(one);
1120  delete one;
1121  delete beta;
1122  return meas;
1123 }
1124 //____________________________________________________________________
1125 TH1* GetDelta(TSeqCollection* l,
1126  const char* which,
1127  Int_t b,
1128  Int_t nEv,
1129  Double_t scale=0,
1130  Double_t scaleE=0,
1131  Double_t lim=1e9)
1132 {
1133  TH1* h = CopyH1(l, Form("b%d_Tr%s_WDist", b, which));
1134  h->Scale(1./nEv);
1135  h->SetYTitle("tracklets/event");
1136  h->SetXTitle("#Delta");
1137 
1138  Bool_t isData = (scale < 1e-6);
1139  Color_t color = (isData ? kGreen+2 : kBlue+2);
1140  TString w(which);
1141  if (w.EqualTo("Data")) SetAttr(h, color, isData ? 21 : 25);
1142  if (w.EqualTo("Inj")) SetAttr(h, color, isData ? 23 : 32);
1143  if (w.EqualTo("Comb")) SetAttr(h, color, 26);
1144  if (w.EqualTo("Prim")) SetAttr(h, kMagenta+2, 28);
1145  if (w.EqualTo("Sec")) SetAttr(h, kMagenta+2, 27);
1146 
1147  if (w.EqualTo("Data")) h->SetTitle("Measured");
1148  if (w.EqualTo("Inj")) h->SetTitle("Injection");
1149  if (w.EqualTo("Comb")) h->SetTitle("Combinatorial");
1150  if (w.EqualTo("Prim")) h->SetTitle("Primaries");
1151  if (w.EqualTo("Sec")) h->SetTitle("Secondaries");
1152 
1153  if (w.EqualTo("Inj")) {
1154  Double_t eR, r = GetDeltaScale(l, b, eR, 1, lim);
1155  Scale(h, r, eR);
1156  h->SetTitle(Form("%5.3f#pm%5.3f #times %s",
1157  r, eR, h->GetTitle()));
1158  }
1159 
1160  if (scale > 1e-6) {
1161  Scale(h, scale, scaleE);
1162  h->SetTitle(Form("%5.3f#pm%5.3f #times %s",
1163  scale, scaleE, h->GetTitle()));
1164  }
1165  return h;
1166 }
1167 
1168 
1169 //____________________________________________________________________
1252 void CorrectOneBin(THStack* stack,
1253  TDirectory* dir,
1254  Int_t b,
1255  TSeqCollection* realList,
1256  TSeqCollection* simList,
1257  Bool_t realComb,
1258  Bool_t simComb,
1259  Bool_t preScale=false,
1260  Bool_t full=false,
1261  Bool_t showOne=false,
1262  UShort_t scaleMode=2,
1263  Double_t lim=1e9)
1264 {
1265  const Double_t k = (realList == simList ? 1 : 1.3);
1266  realList->SetName("realList");
1267  simList ->SetName("simList");
1268  TH2* realMeas = GetH2(realList, Form("b%d_TrData_ZvEtaCutT", b));
1269  TH2* realIPzC = GetH2(realList, "zv");
1270  TH1* realIPz = realIPzC->ProjectionX("realIPz", b+1,b+1, "e");
1271  TH2* simMeas = Coerce(realMeas,
1272  GetH2(simList, Form("b%d_TrData_ZvEtaCutT", b)));
1273  TH2* simIPzC = Coerce(realIPzC, GetH2(simList, "zv"));
1274  TH1* simIPz = simIPzC->ProjectionX("simIPz", b+1,b+1, "e");
1275  TH2* trueGen = Coerce(realMeas,
1276  GetH2(simList, Form("b%d_zvEtaPrimMC", b)));
1277  TH2* trueIPzC = Coerce(realIPzC, GetH2(simList, "zvMCNoPS"));
1278  TH1* trueIPz = simIPzC->ProjectionX("trueIPz", b+1,b+1, "e");
1279  TH2* simBg = Coerce(realMeas, GetBg(simComb, simList, b, 1, lim));
1280  TH2* realBg = Coerce(realMeas, GetBg(realComb,
1281  (realComb ?simList : realList),
1282  b, 1, lim));
1283  Double_t scale = 1;
1284  Double_t scaleE = 0;
1285  TH1* histK = 0;
1286  // PrintH(simBg);
1287  // PrintH(realBg);
1288  if (realComb) {
1289  if (scaleMode == 0) {
1290  // Fixed scale
1291  scale = k;
1292  Scale(realBg, scale, 0);
1293  }
1294  else if (scaleMode == 1) {
1295  // Eta indepedent scale
1296 #if 0
1297  Double_t simRe, realRe;
1298  Double_t realR = GetDeltaScale(realList, b, realRe,
1299  realIPz->GetEntries(), lim);
1300  Double_t simR = GetDeltaScale(simList, b, simRe,
1301  simIPz->GetEntries(), lim);
1302  scale = RatioE(realR, realRe, simR, simRe, scaleE);
1303 #else
1304  scale = GetDeltaScale(realList, simList, b, scaleE,
1305  realIPz->GetEntries(),
1306  simIPz->GetEntries(),
1307  lim);
1308 #endif
1309  Scale(realBg, scale, scaleE);
1310  }
1311  else if (scaleMode == 2) {
1312 #if 0
1313  // Eta dependent scale
1314  TH1* tmpReal = GetDeltaDist(realList, b, realIPz->GetEntries(), lim);
1315  TH1* tmpSim = GetDeltaDist(simList, b, simIPz ->GetEntries(), lim);
1316  tmpReal->Divide(tmpSim);
1317  tmpReal->SetName("k");
1318  tmpReal->SetTitle("Background scale");
1319  tmpReal->SetMarkerColor(tmpReal->GetMarkerColor());
1320  tmpReal->SetMarkerStyle(tmpReal->GetMarkerStyle()+4);
1321  tmpReal->SetLineColor(tmpReal->GetLineColor());
1322  histK = tmpReal;
1323  Scale(realBg, tmpReal);
1324 #else
1325  histK = GetDeltaDist(realList, simList, b,
1326  realIPz->GetEntries(),
1327  simIPz ->GetEntries(), lim);
1328  TF1* kfit = new TF1("kfit", "pol0");
1329  histK->Fit(kfit, "Q0+");
1330  scale = kfit->GetParameter(0);
1331  scaleE = kfit->GetParError(0);
1332  Scale(realBg, histK);
1333 #endif
1334  }
1335  else {
1336  scale = 1;
1337  scaleE = 0;
1338  }
1339  }
1340  if (!CheckConsistency(realMeas, simMeas)) {
1341  Warning("CorrectOneBin", "Real data (%s) and simulated data (%s) done "
1342  "with different binning", realList->GetName(), simList->GetName());
1343  return;
1344  }
1345  realMeas->SetName("realMeas");
1346  simMeas ->SetName("simMeas");
1347  realBg ->SetName("realBg");
1348  simBg ->SetName("simBg");
1349  trueGen ->SetName("trueGen");
1350  realMeas->SetTitle("Measured (real)");
1351  realBg ->SetTitle("Background (real)");
1352  simMeas ->SetTitle("Measured (simulated)");
1353  simBg ->SetTitle("Background (simulated)");
1354  trueGen ->SetTitle("Generated");
1355  realMeas->SetDirectory(0);
1356  simMeas ->SetDirectory(0);
1357  trueGen ->SetDirectory(0);
1358  realIPz ->SetDirectory(0);
1359  simIPz ->SetDirectory(0);
1360  trueIPz ->SetDirectory(0);
1361 
1362  // Scale IPz histograms
1363  Double_t dummy;
1364  Double_t ipzMin = realMeas->GetYaxis()->GetXmin();
1365  Double_t ipzMax = realMeas->GetYaxis()->GetXmax();
1366  Double_t realNev = Integrate(realIPz,ipzMin,ipzMax,dummy);
1367  Double_t simNev = Integrate(simIPz, ipzMin,ipzMax,dummy);
1368  Double_t trueNev = Integrate(trueIPz,ipzMin,ipzMax,dummy);
1369  realIPz->SetMarkerColor(kRed+2); realIPz->SetMarkerStyle(20);
1370  simIPz ->SetMarkerColor(kBlue+2); simIPz ->SetMarkerStyle(21);
1371  trueIPz->SetMarkerColor(kBlack); trueIPz->SetMarkerStyle(24);
1372 #if 0
1373  // This is done in Rubens code
1374  for (Int_t i = 1; i <= realIPz->GetNbinsX(); i++) {
1375  if (simIPz->GetBinContent(i) < 1e-6) {
1376  realIPz->SetBinContent(i,0);
1377  realIPz->SetBinError (i,0);
1378  }
1379  }
1380 #endif
1381 
1382  // Either scale fully here - up front
1383  if (preScale) {
1384  CorrectForIPz(realMeas, realIPz, full);
1385  CorrectForIPz(realBg, (realComb ? simIPz : realIPz), full);
1386  CorrectForIPz(simMeas, simIPz, full);
1387  CorrectForIPz(simBg, simIPz, full);
1388  CorrectForIPz(trueGen, trueIPz, full);
1389  }
1390  // Or first scale by number of events in each sample
1391  else {
1392  // Printf("Number of real events: %f", realNev);
1393  realMeas->Scale(1/realNev);
1394  realBg ->Scale(1/(realComb ? simNev : realNev));
1395  realIPz ->Scale(1/realNev);
1396  simMeas ->Scale(1/simNev);
1397  simBg ->Scale(1/simNev);
1398  simIPz ->Scale(1/simNev);
1399  trueGen ->Scale(1/trueNev);
1400  trueIPz ->Scale(1/trueNev);
1401  }
1402 
1403  // Calculate real signal as measured minus background
1404  TH2* realSig = static_cast<TH2*>(realMeas->Clone("realSig"));
1405  realSig->SetTitle("Signal (real)");
1406  realSig->SetDirectory(0);
1407  CorrectSignal(realComb, realSig, realBg, simMeas);
1408 
1409  // Calculate simulated signal as measured minus background
1410  TH2* simSig = static_cast<TH2*>(simMeas->Clone("simSig"));
1411  simSig->SetTitle("Signal (simulated)");
1412  simSig->SetDirectory(0);
1413  // CorrectSignal(simComb, simSig, simBg, simMeas);
1414  simSig->Add(simBg, -1);
1415 
1416  // Calculate alpha as primary over simulated signal
1417  TH2* alpha = static_cast<TH2*>(trueGen->Clone("alpha"));
1418  alpha->SetTitle("#alpha");
1419  alpha->SetDirectory(0);
1420  alpha->Divide(simSig);
1421  // Create fiducial cut as cut on alpha
1422  TH2* fiducial = static_cast<TH2*>(alpha->Clone("fiducial"));
1423  fiducial->SetTitle("Fiducial cut");
1424  fiducial->SetDirectory(0);
1425  for (Int_t i = 1; i <= fiducial->GetNbinsX(); i++) {
1426  for (Int_t j = 1; j <= fiducial->GetNbinsY(); j++) {
1427  Double_t a = fiducial->GetBinContent(i,j);
1428  fiducial->SetBinError(i,j,0);
1429  fiducial->SetBinContent(i,j, (a > 0 && a <= 2.5));
1430  }
1431  }
1432  alpha ->Multiply(fiducial);
1433 
1434  // Create result as alpha times real signal
1435  TH2* result = static_cast<TH2*>(realSig->Clone("result"));
1436  result->SetTitle("Result");
1437  result->SetDirectory(0);
1438  result->Multiply(alpha);
1439 
1440  // The average mode
1441  UShort_t mode = (preScale ? (full ? 2 : 1) : 0);
1442  // Calculate the primary dN/deta
1443  TH1* truth = Avg(trueGen, trueIPz, mode, "truth");
1444  truth->SetYTitle(dndetaTitle);
1445  SetAttr(truth, cc[b%11], 24, 1.5, 0, 3, 2);
1446 
1447  // Calculate the real dN/deta
1448  TH1* dndeta = Avg(result, realIPz, mode, "dndeta");
1449  dndeta->SetYTitle(dndetaTitle);
1450  SetAttr(dndeta, cc[b%11], 20, 1.2, 0, 1, 1);
1451  if (full) {
1452  for (Int_t i = 1; i <= dndeta->GetNbinsX(); i++) {
1453  dndeta->SetBinError(i,1./3*dndeta->GetBinError(i));
1454  }
1455  }
1456  Double_t lim = 0.5;
1457  TF1* f = new TF1("ff", "pol0", -lim, lim);
1458  TFitResultPtr r = dndeta->Fit(f, "QN0RS", "", -lim, lim);
1459  Printf("dNch/deta %20s: %6.1f +/- %6.1f (%5.2f - %5.3f+/-%5.3f)",
1460  dir->GetName(),
1461  r->Parameter(0), r->ParError(0), r->Chi2()/r->Ndf(),
1462  scale, scaleE);
1463 
1464 
1465  // Multiply histograms by fiducial cut
1466  realMeas->Multiply(fiducial);
1467  simMeas ->Multiply(fiducial);
1468  realSig ->Multiply(fiducial);
1469  simSig ->Multiply(fiducial);
1470  realBg ->Multiply(fiducial);
1471  simBg ->Multiply(fiducial);
1472  realSig ->Multiply(fiducial);
1473  simSig ->Multiply(fiducial);
1474  TH1* realAvgMeas = Avg(realMeas,realIPz, mode, "realAvgMeas",result);
1475  TH1* simAvgMeas = Avg(simMeas, simIPz, mode, "simAvgMeas", result);
1476  TH1* realAvgBg = Avg(realBg, realIPz, mode, "realAvgBg", result);
1477  TH1* simAvgBg = Avg(simBg, simIPz, mode, "simAvgBg", result);
1478  TH1* realAvgSig = Avg(realSig, realIPz, mode, "realAvgSig", result);
1479  TH1* simAvgSig = Avg(simSig, simIPz, mode, "simAvgSig", result);
1480  SetAttr(realAvgMeas, kGreen+2, 21);
1481  SetAttr(realAvgBg, kGreen+2, 22);
1482  SetAttr(realAvgSig, kGreen+2, 23);
1483  SetAttr(simAvgMeas, kRed+2, 22);
1484  SetAttr(simAvgBg, kBlue+2, 21);
1485  SetAttr(simAvgSig, kBlue+2, 22);
1486  SetAttr(histK, kGreen+2, 23);
1487 
1488  THStack* summary = new THStack("summary", dndeta->GetYaxis()->GetTitle());
1489  summary->Add(dndeta, "e2");
1490  summary->Add(truth, "E");
1491  summary->Add(realAvgMeas);
1492  summary->Add(simAvgMeas);
1493  summary->Add(realAvgBg);
1494  summary->Add(simAvgBg);
1495  summary->Add(realAvgSig);
1496  summary->Add(simAvgSig);
1497  if (histK) {
1498  // histK->Scale(10);
1499  // histK->SetTitle(Form("%s\\times10",histK->GetTitle()));
1500  summary->Add(histK);
1501  }
1502 
1503  if (preScale) {
1504  // Scale vertex distributions
1505  realIPz->Scale(1/realNev);
1506  simIPz ->Scale(1/simNev);
1507  trueIPz->Scale(1/trueNev);
1508  }
1509  if (showOne){
1510  TCanvas* c = new TCanvas(Form("c%02d",b),
1511  Form("Centrality bin %d",b),
1512  1200, 1000);
1513  c->Divide(4,3); //,0,0);
1514  c->cd(1); realMeas->DrawCopy("lego2 e");
1515  c->cd(2); realBg ->DrawCopy("lego2 e");
1516  c->cd(3); realSig ->DrawCopy("lego2 e");
1517  c->cd(4); realIPz ->DrawCopy();
1518  c->cd(4); simIPz ->DrawCopy("same");
1519  c->cd(4); trueIPz ->DrawCopy("same");
1520  c->cd(5); simMeas ->DrawCopy("lego2 e");
1521  c->cd(6); simBg ->DrawCopy("lego2 e");
1522  c->cd(7); simSig ->DrawCopy("lego2 e");
1523  c->cd(8); alpha ->DrawCopy("lego2 e");
1524  c->cd(9); fiducial->DrawCopy("lego2 e");
1525  c->cd(10); trueGen ->DrawCopy("lego2 e");
1526  c->cd(11); result ->DrawCopy("lego2 e");
1527  c->cd(12); summary ->DrawClone("nostack");
1528  c->GetPad(4)->Modified();
1529  c->GetPad(12)->Modified();
1530  c->GetPad(4)->BuildLegend();
1531  c->GetPad(12)->BuildLegend();
1532  TLatex* ltx = new TLatex(.52, .3,
1533  Form("d#it{N}_{ch}/d#eta|_{|#eta|<%3.1f}"
1534  "=%6.1f#pm%6.1f (#chi^{2}/#nu=%5.2f)",
1535  lim, r->Parameter(0),
1536  r->ParError(0), r->Chi2()/r->Ndf()));
1537  ltx->SetTextAlign(23);
1538  ltx->SetNDC(true);
1539  ltx->SetTextFont(42);
1540  ltx->SetTextSize(0.035);
1541  c->cd(12); ltx->Draw();
1542  ltx->DrawLatex(.52,.25,
1543  Form("%sk%s=%5.3f#pm%5.3f",
1544  histK ? "#LT" : "",
1545  histK ? "#GT" : "",
1546  scale,
1547  scaleE));
1548  }
1549  stack->Add(truth);
1550  stack->Add(dndeta, "e2");
1551 
1552  if (!dir) return;
1553 
1554  dir->cd();
1555  result ->Write();
1556  summary->Write();
1557  dndeta ->Write();
1558  if (histK) histK ->Write();
1559  (new TParameter<double>("scale", scale)) ->Write();
1560  (new TParameter<double>("scaleE", scaleE))->Write();
1561 
1562  TDirectory* det = dir->mkdir("details");
1563  det->cd();
1564  if (!preScale) {
1565  CorrectForIPz(realMeas, realIPz, full);
1566  CorrectForIPz(realBg, (realComb ? simIPz : realIPz), full);
1567  CorrectForIPz(simMeas, simIPz, full);
1568  CorrectForIPz(simBg, simIPz, full);
1569  CorrectForIPz(trueGen, trueIPz, full);
1570  }
1571  realMeas ->Write();
1572  realBg ->Write();
1573  realSig ->Write();
1574  simMeas ->Write();
1575  simBg ->Write();
1576  simSig ->Write();
1577  trueGen ->Write();
1578  alpha ->Write();
1579  fiducial ->Write();
1580  if (histK) histK->Write();
1581  THStack* deltas = new THStack("deltas", "#Delta distributions");
1582  deltas->Add(GetDelta(realList,"Data",b,realNev,0, 0, lim));
1583  deltas->Add(GetDelta(simList, "Data",b,simNev, scale,scaleE,lim));
1584  deltas->Add(GetDelta(realList,"Inj", b,realNev,0, 0, lim));
1585  deltas->Add(GetDelta(simList, "Inj", b,simNev, scale,scaleE,lim));
1586  deltas->Add(GetDelta(simList, "Comb",b,simNev, scale,scaleE,lim));
1587  deltas->Add(GetDelta(simList, "Prim",b,trueNev,scale,scaleE,lim));
1588  deltas->Add(GetDelta(simList, "Sec", b,trueNev,scale,scaleE,lim));
1589  TH1* sum = GetDelta(simList, "Prim", b,trueNev,0, 0, lim);
1590  sum->Add(GetDelta(simList, "Sec", b,trueNev,0, 0, lim));
1591  sum->Add(GetDelta(simList, "Comb", b,simNev, 0, 0, lim));
1592  sum->SetTitle(Form("%5.3f#pm%5.3f #times (Prim+Sec+Comb)",scale,scaleE));
1593  Scale(sum,scale,scaleE);
1594  SetAttr(sum, kPink+2, 30);
1595  deltas->Add(sum);
1596  deltas->Write();
1597 
1598  TDirectory* avs = dir->mkdir("averages");
1599  avs->cd();
1600  realAvgMeas->Write();
1601  realAvgSig ->Write();
1602  realAvgBg ->Write();
1603  simAvgMeas ->Write();
1604  simAvgSig ->Write();
1605  simAvgBg ->Write();
1606  realIPz ->Write();
1607  simIPz ->Write();
1608  trueIPz ->Write();
1609  truth ->Write();
1610 
1611 }
1612 /* @} */
1613 //====================================================================
1625 void
1626 ShowStats(TList* list, const char* title, const char* path)
1627 {
1628  Printf("=== %-15s - %s ===", title, path);
1629  TH1* tmp = GetH1(list, "hStat");
1630  if (!tmp) {
1631  Warning("", "No stats");
1632  return;
1633  }
1634  TH1* stats = static_cast<TH1*>(tmp->Clone());
1635  stats->SetDirectory(0);
1636  TAxis* axis = stats->GetXaxis();
1637  for (Int_t i = 1; i <= 4; i++) {
1638  Printf("%31s: %f", axis->GetBinLabel(i), stats->GetBinContent(i));
1639  }
1640  Double_t scale = stats->GetBinContent(4);
1641  stats->Scale(1/scale);
1642  for (Int_t i = 6; i <= 23; i++) {
1643  Printf("%31s: %f", axis->GetBinLabel(i), stats->GetBinContent(i));
1644  }
1645 }
1646 /* @} */
1647 
1648 //____________________________________________________________________
1671  UShort_t scaleMode,
1672  const char* realFileName,
1673  const char* simFileName,
1674  Int_t nBins=9,
1675  const char* otherFileName="")
1676 {
1677  TFile* realFile = TFile::Open(realFileName, "READ");
1678  TFile* simFile = TFile::Open(simFileName, "READ");
1679  TFile* otherFile = 0;
1680  if (otherFileName && otherFileName[0] != '\0')
1681  otherFile = TFile::Open(otherFileName, "READ");
1682  TList* realList = static_cast<TList*>(realFile->Get("clist"));
1683  TList* simList = static_cast<TList*>(simFile ->Get("clist"));
1684  TObjArray* otherList = 0;
1685  if (otherFile)
1686  otherList =static_cast<TObjArray*>(otherFile->Get("TObjArray"));
1687 
1688  ShowStats(realList, "Real data", realFile->GetName());
1689  ShowStats(simList, "Simulated data",simFile ->GetName());
1690 
1691  TH1* realCent = GetH1(realList, "EvCentr");
1692  TH1* simCent = GetH1(simList, "EvCentr");
1693  if (!CheckConsistency(realCent, simCent)) {
1694  Warning("SimpleCorrection","Inconsistent centralities");
1695  return;
1696  }
1697  TH2* realIpz = GetH2(realList, "b0_TrData_ZvEtaCutT");
1698  TH2* simIpz = GetH2(simList, "b0_TrData_ZvEtaCutT");
1699  if (!CheckConsistency(realIpz, simIpz)) {
1700  Warning("SimpleCorrection", "Inconsistent IPz or eta axis");
1701  return;
1702  }
1703 
1704 
1705  Double_t deltaLim = 1e9;
1706  Bool_t realComb = flags & 0x1;
1707  Bool_t simComb = flags & 0x2;
1708  Bool_t preScale = flags & 0x4;
1709  Bool_t full = flags & 0x8;
1710  Bool_t showOne = flags & 0x20;
1711  // UShort_t scaleMode = (flags & 0x40 ? 0 : (flags & 0x80 ? 1 : 2));
1712  Printf("Real background: %s\n"
1713  "Simulated background: %s\n"
1714  "Pre-scale: %s\n"
1715  "Full pre-scale: %s\n"
1716  "Show each bin: %s\n"
1717  "Delta upper limit: %g\n"
1718  "Scaling mode: %s",
1719  (realComb ? "comb" : "inj"),
1720  (simComb ? "comb" : "inj"),
1721  (preScale ? "yes" : "no"),
1722  (full ? "yes" : "no"),
1723  (showOne ? "yes" : "no"),
1724  deltaLim,
1725  (scaleMode==0 ?"fixed" :
1726  scaleMode==1 ? "c" :
1727  scaleMode==2 ? "eta-c" :
1728  "unit"));
1729 
1730 
1731  if (flags & 0x10) realList = simList;
1732  TString tit;
1733  tit.Form("Real BG: %s Simulated BG: %s%s",
1734  realComb ? "MC-labels" : "Injection",
1735  simComb ? "MC-labels" : "Injection",
1736  preScale ? " (pre-scaled)" : "");
1737 
1738  TFile* output = TFile::Open(Form("result_0x%x.root", flags&0x3),
1739  "RECREATE");
1740  realCent->Write("cent");
1741  THStack* res = new THStack("result", tit);
1742  TH1* cent = GetH1(realList, "EvCentr");
1743  for (Int_t b = 0; b < nBins; b++) {
1744  Double_t c1 = cent->GetXaxis()->GetBinLowEdge(b+1);
1745  Double_t c2 = cent->GetXaxis()->GetBinUpEdge (b+1);
1746  TString name;
1747  name.Form("cent%03dd%02d_%03dd%02d",
1748  Int_t(c1), Int_t(c1*100)%100,
1749  Int_t(c2), Int_t(c2*100)%100);
1750  TDirectory* dir = output->mkdir(name);
1751  CorrectOneBin(res,dir,b,realList,simList,realComb,simComb,preScale,full,
1752  showOne, scaleMode, deltaLim);
1753 
1754  if (otherList) {
1755  TH1* od = GetH1(otherList, Form("bin%d_DataCorrSignal_Bla",b));
1756  TH1* ot = GetH1(otherList, Form("bin%d_MCTruth",b));
1757  TH1* td = static_cast<TH1*>(od->Clone("otherdNdeta"));
1758  TH1* tt = static_cast<TH1*>(ot->Clone("otherTruth"));
1759  td->SetDirectory(0);
1760  tt->SetDirectory(0);
1761  SetAttr(td, cc[b%11], 21, 1.2, 3004, 1, 1);
1762  SetAttr(tt, cc[b%11], 25, 1.4, 0, 7, 3);
1763  if (td->GetMinimum() < 1e-6) { td->SetMinimum(); td->SetMaximum(); }
1764  if (tt->GetMinimum() < 1e-6) { tt->SetMinimum(); tt->SetMaximum(); }
1765  res->Add(td, "e2");
1766  res->Add(tt);
1767  dir->cd();
1768  td->Write();
1769  tt->Write();
1770  }
1771  }
1772  // TCanvas* all = new TCanvas("all", "all");
1773  // res->Draw("nostack");
1774  res->SetMaximum(1.2*res->GetMaximum("nostack"));
1775  output->cd();
1776  res->Write();
1777  output->Write();
1778 
1779 }
1780 /* @} */
1781 //
1782 // EOF
1783 //
Bool_t CheckAxisNBins(const char *which, const TAxis *a1, const TAxis *a2)
Int_t color[]
const char * dndetaTitle
Definition: SimpleCorrect.C:31
Bool_t CheckAxisLabels(const char *which, const TAxis *a1, const TAxis *a2)
double Double_t
Definition: External.C:58
const Bool_t kSimpleCorrectLoaded
Definition: SimpleCorrect.C:26
const char * title
Definition: MakeQAPdf.C:26
Bool_t CheckAxisBins(const char *which, const TAxis *a1, const TAxis *a2)
TH1 * Avg2(TH2 *h, TH1 *ipz, UShort_t mode, const char *name)
TList * list
TH1 * GetH1(TSeqCollection *l, const char *name)
TCanvas * c
Definition: TestFitELoss.C:172
TH1 * SetAttr(TH1 *h, Color_t color, Style_t marker=20, Double_t size=1., Style_t fill=0, Style_t line=1, Width_t width=1)
void ShowStats(TList *list, const char *title, const char *path)
Double_t Integrate(TH1 *h, Double_t min, Double_t max, Double_t &err)
TH2 * GetH2(TSeqCollection *l, const char *name)
TH2 * CorrectSignal(Bool_t comb, TH2 *meas, TH2 *bg, TH2 *simMeas)
Bool_t CheckConsistency(const TH1 *h1, const TH1 *h2)
TH1 * GetDeltaDist(TH2 *deltaReal, TH2 *deltaSim, Double_t nEvReal, Double_t nEvSim, Double_t lim=1e9)
int Int_t
Definition: External.C:63
TH1 * Avg(TH2 *h, TH1 *ipz, UShort_t mode, const char *name, TH2 *other=0)
TH2 * CopyH2(TSeqCollection *l, const char *name, const char *newName=0)
Double_t RatioE(Double_t n, Double_t en, Double_t d, Double_t ed, Double_t &er)
TH2 * GetBg(Bool_t comb, TSeqCollection *list, Int_t b, Double_t nEv=1, Double_t lim=1e9)
TObject * CopyO(TSeqCollection *l, const char *name, const char *newName, TClass *cls)
Definition: SimpleCorrect.C:83
Int_t mode
Definition: anaM.C:40
Bool_t CheckAxisLimits(const char *which, const TAxis *a1, const TAxis *a2)
void PrintH(TH2 *h, Int_t prec=2)
void SimpleCorrect(UShort_t flags, UShort_t scaleMode, const char *realFileName, const char *simFileName, Int_t nBins=9, const char *otherFileName="")
Definition: External.C:220
TH2 * Scale(TH2 *h, TH1 *g)
unsigned short UShort_t
Definition: External.C:28
TObject * GetO(TSeqCollection *l, const char *name, TClass *cls)
Definition: SimpleCorrect.C:61
TH1 * CopyH1(TSeqCollection *l, const char *name, const char *newName=0)
Bool_t CheckAxis(const char *which, const TAxis *a1, const TAxis *a2, Bool_t alsoLbls)
Double_t GetDeltaScale(TH1 *deltaRec, TH1 *deltaBg, Double_t &eR, Double_t nEvRec=1, Double_t nEvBg=1, Double_t lim=1e9)
bool Bool_t
Definition: External.C:53
TH1 * Coerce(const TH1 *templ, TH1 *src)
TH2 * CorrectForIPz(TH2 *in, TH1 *ipz, Bool_t full=true)
Definition: External.C:196
void CorrectOneBin(THStack *stack, TDirectory *dir, Int_t b, TSeqCollection *realList, TSeqCollection *simList, Bool_t realComb, Bool_t simComb, Bool_t preScale=false, Bool_t full=false, Bool_t showOne=false, UShort_t scaleMode=2, Double_t lim=1e9)
const Color_t cc[]
Definition: SimpleCorrect.C:35
TH1 * GetDelta(TSeqCollection *l, const char *which, Int_t b, Int_t nEv, Double_t scale=0, Double_t scaleE=0, Double_t lim=1e9)