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