AliPhysics  cda3415 (cda3415)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliTrackletWeights.C
Go to the documentation of this file.
1 
11 #ifndef ALITRACKLETWEIGHTS_C
12 #define ALITRACKLETWEIGHTS_C
13 #include <TNamed.h>
14 #include <map>
15 #ifndef __CINT__
16 # include <TH1.h>
17 # include <TH2.h>
18 # include <TH3.h>
19 # include <TAxis.h>
20 # include <TList.h>
21 # include <TBrowser.h>
22 # include <TCanvas.h>
23 # include <THStack.h>
24 # include <TROOT.h>
25 # include <TClass.h>
26 # include "AliAODTracklet.C"
27 # include <TParameter.h>
28 #else
29 class TH1D;
30 class TH2D;
31 class TH1;
32 class TH2;
33 class TH3;
34 class TAxis;
35 class TList;
36 class TBrowser;
37 class THStack;
38 class AliAODTracklet;
39 class TCanvas; // Autoload
40 class TVirtualPad;
41 #endif
42 
43 //====================================================================
44 #ifndef __CINT__
45 namespace {
46  template <typename T>
47  T GetPar(const char* name, TCollection* l)
48  {
49  TObject* o = l->FindObject(name);
50  if (!o) {
51  Warning("GetPar", "Didn't find parameter %s in %s",
52  name, l->GetName());
53  return T();
54  }
55  TClass* cls = TParameter<T>::Class();
56  if (!o->IsA()->InheritsFrom(cls)) {
57  Warning("GetPar", "Object %s is a %s, not a %s",
58  name, o->ClassName(), cls->GetName());
59  return T();
60  }
61  TParameter<T>* p = static_cast<TParameter<T>*>(o);
62  return p->GetVal();
63  }
64 }
65 #endif
66 
67 //====================================================================
72 {
73 public:
77  enum ECalc {
120  };
126  : TNamed(),
127  fCalc(kProduct),
128  fMask(0xFF),
129  fVeto(0x0)
130  {}
137  AliTrackletBaseWeights(const char* name,
138  const char* title="Sim. tracklet weights")
139  : TNamed(name,title),
140  fCalc(kProduct),
141  fMask(0xFF),
142  fVeto(0x0),
143  fInverse(false)
144  {}
151  : TNamed(o),
152  fCalc (o.fCalc),
153  fMask (o.fMask),
154  fVeto (o.fVeto),
155  fInverse(o.fInverse)
156  {}
161 
170  {
171  if (&o == this) return *this;
172  TNamed::operator=(o);
173  fCalc = o.fCalc;
174  fMask = o.fMask;
175  fVeto = o.fVeto;
176  fInverse = o.fInverse;
177  return *this;
178  }
226  void SetCalc(UChar_t mode=kProduct) { fCalc = mode; }
232  void SetMask(UChar_t mask) { fMask = mask; }
238  void SetVeto(UChar_t veto) { fVeto = veto; }
246  void SetInverse(Bool_t inv) { fInverse = inv; }
254  Bool_t CheckTracklet(const AliAODTracklet* tracklet) const
255  {
256  UChar_t flags = tracklet->GetFlags();
257  if (fMask != 0xFF && (fMask & flags) == 0) {
258  // Info("LookupWeight", "Tracklet 0x%02x does not fullfill mask 0x%02x",
259  // flags, fMask);
260  return false;
261  }
262  if ((fVeto & flags) != 0) {
263  // Info("LookupWeight", "Tracklet 0x%02x vetoed by 0x%02x",flags, fVeto);
264  return false;
265  }
266  return true;
267  }
279  Double_t cent,
280  Double_t ipz,
281  TH2* corr=0) const
282  {
283  if (!CheckTracklet(tracklet)) return 1;
284  Double_t w = CalcWeight(tracklet, cent, ipz, corr);
285  if (fInverse) w = 1/w;
286  return w;
287  }
299  virtual Double_t CalcWeight(AliAODTracklet* tracklet,
300  Double_t cent,
301  Double_t ipz,
302  TH2* corr) const = 0;
308  virtual TCollection* Store(TCollection* parent)
309  {
310  TList* top = new TList;
311  top->SetName(GetName());
312  top->SetOwner(true);
313  parent->Add(top);
314  top->Add(new TParameter<int>("mask", fMask, 'f'));
315  top->Add(new TParameter<int>("veto", fVeto, 'f'));
316  top->Add(new TParameter<int>("calc", fCalc, 'f'));
317  return top;
318  }
327  {
328  TCollection* top = static_cast<TCollection*>(in->FindObject(GetName()));
329  if (!top) {
330  Warning("Retrieve",
331  "Collection %s not found in %s", GetName(), in->GetName());
332  in->ls();
333  return 0;
334  }
335  fCalc = GetPar<int>("calc", top);
336  fMask = GetPar<int>("mask", top);
337  fVeto = GetPar<int>("veto", top);
338  return top;
339  }
343  Bool_t IsFolder() const { return true; }
349  virtual void Print(Option_t* option="") const //*MENU*
350  {
351  gROOT->IndentLevel();
352  Printf("%s : %s", ClassName(), GetName());
353  Printf(" Weight calculation: %s",
354  fCalc == kProduct ? "product" :
355  fCalc == kSquare ? "square" :
356  fCalc == kSum ? "sum" : "average");
357  Printf(" Tracklet mask: 0x%02x", fMask);
358  Printf(" Tracklet veto: 0x%02x", fVeto);
359  }
360  UChar_t fCalc; // Whether the square of the weight is calculated
361  UChar_t fMask; // Which partiles to take
362  UChar_t fVeto; // Which particles not to take
363  Bool_t fInverse; // If true, do 1/w
364  ClassDef(AliTrackletBaseWeights,1); // Base class of weights
365 };
366 
367 //====================================================================
374 {
375 public:
379  enum EMode {
381  kUp = (0x1) << 14,
383  kDown = (0x2) << 14,
385  kDisabled = (0x4) << 14
386  };
387 
389  typedef std::map<short,TH1D*> PdgMap;
396  fPt(0),
397  fAbundance(),
398  fStrangeness()
399  {}
406  AliTrackletPtPidStrWeights(const char* name,
407  const char* title="Sim. tracklet weights");
426 
437  virtual Double_t CalcWeight(AliAODTracklet* tracklet,
438  Double_t cent,
439  Double_t ipz,
440  TH2* corr=0) const;
451  {
452  AddPdgWeight(fAbundance, pdg, h, mode);
453  }
464  {
465  AddPdgWeight(fStrangeness, pdg, h, mode);
466  }
475  Bool_t SetPtWeight(const TH2D* h, UShort_t mode=0);
489  {
490  SetPdgMode(fAbundance, pdg, mode);
491  }
504  {
505  SetPdgMode(fStrangeness, pdg, mode);
506  }
517  void SetPtMode(UShort_t mode);
523  void Draw(Option_t* option=""); //*MENU*
529  void Print(Option_t* option="") const; //*MENU*
556  {
557  return GetPdgWeight(fAbundance, apdg, cent);
558  }
569  {
570  return GetPdgWeight(fStrangeness, apdg, cent);
571  }
572 
573 protected:
585  virtual Double_t CalcWeight(Double_t pT, Short_t pdg, Double_t cent) const;
594  TH1D* GetPdgHist(const PdgMap& m, Short_t pdg) const;
604  virtual Double_t GetPdgWeight(const PdgMap& m,
605  UShort_t apdg,
606  Double_t cent) const;
617  Bool_t AddPdgWeight(PdgMap& map, Short_t pdg, const TH1D* w, UShort_t mode=0);
629  void SetPdgMode(PdgMap& map, Short_t pdg, UShort_t mode);
637  void StoreMap(TCollection* parent, const char* name, PdgMap& m);
647  Bool_t RetrieveMap(TCollection* parent, const char* name, PdgMap& m);
653  void ModStack(THStack* stack);
661  void PrintMap(const PdgMap& m, const char* name, Option_t* options="") const;
668  void PrintHist(const char* tag, TH1* h) const;
669 
670  TH2D* fPt; // Weight by centrality and pT
671  PdgMap fAbundance; // Map for abundance weight
672  PdgMap fStrangeness; // Map for strangeness weight
673  ClassDef(AliTrackletPtPidStrWeights,3); // Weighs for tracklet analysis
674 };
675 
676 //____________________________________________________________________
678  const char* title)
679  : AliTrackletBaseWeights(name, title),
680  fPt(0),
681  fAbundance(),
682  fStrangeness()
683 {}
684 
685 //____________________________________________________________________
689  fPt(0),
690  fAbundance(),
691  fStrangeness()
692 {
693  UInt_t mask = kUp|kDown|kDisabled;
694  SetPtWeight(o.fPt, o.fPt->TestBits(mask));
695  for (PdgMap::const_iterator i = o.fAbundance.begin();
696  i != o.fAbundance.end(); ++i)
697  AddAbundanceWeight(i->first, i->second, i->second->TestBits(mask));
698  for (PdgMap::const_iterator i = o.fStrangeness.begin();
699  i != o.fStrangeness.end(); ++i)
700  AddStrangenessWeight(i->first, i->second, i->second->TestBits(mask));
701 }
702 //____________________________________________________________________
705 {
706  if (&o == this) return *this;
708  fName = o.fName;
709  fTitle = o.fTitle;
710  UInt_t mask = kUp|kDown|kDisabled;
711  SetPtWeight(o.fPt, o.fPt->TestBits(mask));
712  for (PdgMap::const_iterator i = o.fAbundance.begin();
713  i != o.fAbundance.end(); ++i)
714  AddAbundanceWeight(i->first, i->second, i->second->TestBits(mask));
715  for (PdgMap::const_iterator i = o.fStrangeness.begin();
716  i != o.fStrangeness.end(); ++i)
717  AddStrangenessWeight(i->first, i->second, i->second->TestBits(mask));
718  return *this;
719 }
720 //____________________________________________________________________
722  Short_t pdg,
723  const TH1D* w,
724  UShort_t mode)
725 {
726  UShort_t apdg = TMath::Abs(pdg);
727  TH1D* copy = static_cast<TH1D*>(w->Clone(Form("w%d", apdg)));
728  copy->SetDirectory(0);
729  copy->SetXTitle("Centrality [%]");
730  copy->SetBit(mode);
731  PdgMap::iterator i = m.find(apdg);
732  if (i != m.end()) {
733  Warning("AddPdgWeight", "Replacing weight for %d", pdg);
734  delete i->second;
735  i->second = 0;
736  m.erase(i);
737  }
738  m[apdg] = copy;
739  return true;
740 }
741 
742 //____________________________________________________________________
744 {
745  if (fPt) {
746  delete fPt;
747  fPt = 0;
748  }
749  if (!h) return false;
750  fPt = static_cast<TH2D*>(h->Clone("centPt"));
751  fPt->SetDirectory(0);
752  fPt->SetXTitle("Centrality [%]");
753  fPt->SetYTitle("#it{p}_{T} [GeV]");
754  fPt->SetBit(mode);
755  return true;
756 }
757 //____________________________________________________________________
759 {
760  if (fPt) fPt->SetBit(mode);
761 }
762 //____________________________________________________________________
764  Short_t pdg,
765  UShort_t mode)
766 {
767  TH1* h = GetPdgHist(m, pdg);
768  if (h) h->SetBit(mode);
769 }
770 //____________________________________________________________________
772 {
773  UShort_t apdg = TMath::Abs(pdg);
774  PdgMap::const_iterator i = m.find(apdg);
775  if (i == m.end()) return 0;
776  return i->second;
777 }
778 
779 //____________________________________________________________________
781  UShort_t apdg,
782  Double_t cent) const
783 {
784  if (m.size() < 1) return 1;
785  TH1D* h = GetPdgHist(m, apdg);
786  if (!h || h->TestBit(kDisabled)) return 1;
787 
788 #if 1
789  Int_t bC = h->GetXaxis()->FindBin(cent);
790  if (bC < 1 || bC > h->GetNbinsX()) return 1;
791 #else
792  Int_t bC = 1;
793 #endif
794  Double_t add = (h->TestBit(kUp) ? +1 :
795  h->TestBit(kDown) ? -1 : 0) * h->GetBinError(bC);
796  return h->GetBinContent(bC) + add;
797 }
798 //____________________________________________________________________
799 Double_t
801  Short_t pdg,
802  Double_t cent) const
803 {
804  Double_t w = 1;
805  if (fPt && !fPt->TestBit(kDisabled)) {
806  Int_t bC = fPt->GetXaxis()->FindBin(cent);
807  Int_t bpT = fPt->GetYaxis()->FindBin(pT);
808  if (bC >= 1 && bC <= fPt->GetXaxis()->GetNbins() &&
809  bpT >= 1 && bpT <= fPt->GetYaxis()->GetNbins()) {
810  Double_t fac = (pT >= 0.05 ? 1 :
811  (fPt->TestBit(kUp) ? 1.3 :
812  fPt->TestBit(kDown) ? 0.7 : 1));
813  w *= fac*fPt->GetBinContent(bC, bpT);
814  }
815  }
816  UShort_t apdg = TMath::Abs(pdg);
817 
818  w *= GetPdgWeight(fAbundance, apdg, cent);
819  w *= GetPdgWeight(fStrangeness, apdg, cent);
820  // Printf("Weight of pT=%6.3f pdg=%5d cent=%5.1f -> %f", pT, pdg, cent, w);
821  return w;
822 }
823 
824 //____________________________________________________________________
825 Double_t
827  Double_t cent,
828  Double_t ipz,
829  TH2* corr) const
830 {
831 #if 0
832  if (!tracklet->IsSimulated()) {
833  Warning("LookupWeight", "Not a simulated tracklet");
834  return 1;
835  }
836 #endif
837  Double_t w1 = 1, w2 = 1;
838 
839  // AliAODMCTracklet* mc = static_cast<AliAODMCTracklet*>(tracklet);
840  AliAODTracklet* mc = tracklet;
841  Short_t pdg1 = mc->GetParentPdg();
842  Short_t pdg2 = mc->GetParentPdg(true);
843  if (pdg1>0) w1 = CalcWeight(mc->GetParentPt(), pdg1,cent);
844  if (pdg2>0) w2 = CalcWeight(mc->GetParentPt(true),pdg2,cent);
845  else if (fCalc!=kProduct) w2 = w1;
846 
847  if (corr && mc->IsMeasured()) corr->Fill(w1, w2);
848 
849  switch (fCalc) {
850  case kProduct: return w1 * w2;
851  case kSquare: return TMath::Sqrt(w1 * w2);
852  case kSum: return 1+(w1-1)+(w2-1);
853  case kAverage: return 1+((w1-1)+(w2-1))/2;
854  }
855  return 1;
856 }
857 
858 //____________________________________________________________________
860 {
861  if (!stack || !stack->GetHists()) return;
862  TIter next(stack->GetHists());
863  TH1* hist = 0;
864  Color_t colors[] = { kPink+2, kRed+2, kOrange+2, kYellow+2,
865  kSpring+2, kGreen+2, kTeal+2, kCyan+2,
866  kBlue+2, kViolet+2 };
867  Int_t idx = 0;
868  while ((hist = static_cast<TH1*>(next()))) {
869  Color_t c = colors[idx % 10];
870  idx++;
871  hist->SetLineColor(c);
872  hist->SetMarkerColor(c);
873  hist->SetMarkerStyle(20+idx%20);
874  }
875 }
876 
877 //____________________________________________________________________
880 {
882  if (!top) return 0;
883 
884  if (fPt) {
885  TH2* copy = static_cast<TH2*>(fPt->Clone("centPt"));
886  copy->SetDirectory(0);
887  copy->SetBinContent(0,0,1); // For counting merges
888  top->Add(copy);
889  }
890  StoreMap(top, "abundance", fAbundance);
891  StoreMap(top, "strangeness", fStrangeness);
892  return top;
893 }
894 
895 //____________________________________________________________________
896 void
898  const char* name,
899  PdgMap& m)
900 {
901  TList* top = new TList;
902  top->SetName(name);
903  top->SetOwner(true);
904  parent->Add(top);
905  for (PdgMap::const_iterator i = m.begin(); i != m.end(); ++i) {
906  TH1* copy = static_cast<TH1*>(i->second->Clone(Form("w%d",i->first)));
907  copy->SetDirectory(0);
908  copy->SetBinContent(0,1); // For counting merges
909  top->Add(copy);
910  }
911 }
912 
913 
914 //____________________________________________________________________
917 {
919  if (!top) return 0;
920 
921  fPt = static_cast<TH2D*>(top->FindObject("centPt"));
922  if (!fPt)
923  Warning("Retrieve","centPt histogram not found in %s", GetName());
924  else {
925  fPt->SetDirectory(0);
926  Double_t scale = fPt->GetBinContent(0,0);
927  fPt->Scale(1/scale); // Counting merges
928  // fPt->SetBinContent(0,0,1); // Zero merger count
929  }
930  if (!RetrieveMap(top, "abundance", fAbundance)) return 0;
931  if (!RetrieveMap(top, "strangeness", fStrangeness)) return 0;
932 
933  return top;
934 }
935 //____________________________________________________________________
936 Bool_t
938  const char* name,
939  PdgMap& m)
940 {
941  m.clear();
942  TList* top = static_cast<TList*>(parent->FindObject(name));
943  if (!top) {
944  Warning("RetrieveMap",
945  "Collection %s not found in %s", name, parent->GetName());
946  return true;
947  }
948  TIter next(top);
949  TObject* o = 0;
950  while ((o = next())) {
951  if (!o->IsA()->InheritsFrom(TH1D::Class())) continue;
952  TH1D* copy = static_cast<TH1D*>(o);
953  copy->SetDirectory(0);
954  Double_t scale = copy->GetBinContent(0);
955  copy->Scale(1/scale); // Counting merges
956  // copy->SetBinContent(0,1); // Zero merger count
957  TString nme(copy->GetName());
958  nme.Remove(0,1);
959  Int_t apdg = nme.Atoi();
960  AddPdgWeight(m, apdg, copy, copy->TestBits(kUp|kDown|kDisabled));
961  }
962  return true;
963 }
964 
965 //____________________________________________________________________
966 void
968 {
969  TVirtualPad* master = TVirtualPad::Pad();
970  if (!master) {
971  Warning("Draw", "No current pad to draw in");
972  return;
973  }
974 
975  Int_t nPad = 3;
976  Int_t iPad = 0;
977  if (fAbundance.size() <= 0) nPad--;
978  if (fStrangeness.size() <= 0) nPad--;
979  TString opt(option); opt.ToUpper();
980 
981  master->Divide(nPad,1);
982  if (fPt) {
983  THStack* stack = new THStack(fPt, (opt.Contains("C") ? "x" : "y"));
984  ModStack(stack);
985  stack->SetTitle("#it{p}_{T} weights");
986  master->cd(++iPad);
987  stack->Draw("nostack");
988  stack->GetHistogram()->SetXTitle((opt.Contains("C") ?
989  "Centrality [%]" : "#it{p}_{T}"));
990  master->GetPad(iPad)->BuildLegend();
991  master->GetPad(iPad)->Modified();
992  }
993 
994  if (fAbundance.size() > 0) {
995  THStack* ha = new THStack("abundance", "Abundance weights");
996  for (PdgMap::const_iterator i = fAbundance.begin();
997  i != fAbundance.end(); ++i) {
998  ha->Add(i->second);
999  }
1000  ModStack(ha);
1001  master->cd(++iPad);
1002  ha->Draw("nostack");
1003  ha->GetHistogram()->SetXTitle("Centrality [%]");
1004  master->GetPad(iPad)->BuildLegend();
1005  }
1006 
1007  if (fStrangeness.size() > 0) {
1008  THStack* hs = new THStack("strangeness", "Strangeness weights");
1009  for (PdgMap::const_iterator i = fStrangeness.begin();
1010  i != fStrangeness.end(); ++i) {
1011  hs->Add(i->second);
1012  }
1013  ModStack(hs);
1014  master->cd(++iPad);
1015  hs->Draw("nostack");
1016  hs->GetHistogram()->SetXTitle("Centrality [%]");
1017  master->GetPad(iPad)->BuildLegend();
1018  }
1019  master->Modified();
1020 }
1021 
1022 //____________________________________________________________________
1023 void
1025 {
1027  gROOT->IncreaseDirLevel();
1028  PrintHist("pT", fPt);
1029  PrintMap(fAbundance, "Abundance", option);
1030  PrintMap(fStrangeness, "Strangeness", option);
1031 
1032  gROOT->DecreaseDirLevel();
1033 }
1034 //____________________________________________________________________
1035 void
1036 AliTrackletPtPidStrWeights::PrintHist(const char* tag, TH1* h) const
1037 {
1038  if (!h) return;
1039  gROOT->IndentLevel();
1040  Printf("%10s (%c): %p %s/%s", tag,
1041  h->TestBit(kDisabled) ? '0' :
1042  h->TestBit(kUp) ? '+' :
1043  h->TestBit(kDown) ? '-' : '=',
1044  h, h->GetName(), h->GetTitle());
1045 
1046 }
1047 //____________________________________________________________________
1048 void
1050  const char* name,
1051  Option_t* option) const
1052 {
1053  gROOT->IndentLevel();
1054  Printf("Map of PDG codes: %s", name);
1055  gROOT->IncreaseDirLevel();
1056  for (PdgMap::const_iterator i = m.begin(); i != m.end(); ++i)
1057  PrintHist(Form("%10d", i->first), i->second);
1058 
1059  gROOT->DecreaseDirLevel();
1060 }
1061 
1062 
1063 //====================================================================
1068 {
1069 public:
1076  fHistos(),
1077  fCentAxis()
1078  {}
1085  AliTrackletDeltaWeights(const char* name,
1086  const char* title="Sim. tracklet weights");
1097 
1111  void SetCentAxis(const TAxis& axis);
1118  void SetCentAxis(Int_t n, const Double_t* bins);
1126  void SetCentAxis(Int_t n, Double_t low, Double_t high);
1135  Bool_t SetHisto(Int_t bin, TH1* h);
1143  TH1* FindHisto(Double_t cent) const;
1154  static Int_t FindBin(Double_t value, const TAxis* axis);
1165  virtual Double_t CalcWeight(AliAODTracklet* tracklet,
1166  Double_t cent,
1167  Double_t ipZ,
1168  TH2* corr=0) const;
1169 
1179  TH1* Project(TH1* h, char which, const char* name);
1185  void Draw(Option_t* option=""); //*MENU*
1191  void Print(Option_t* option="") const; //*MENU*
1202  void DrawStack(TVirtualPad* pad,
1203  Int_t nPad,
1204  THStack* stack,
1205  const char* xtitle,
1206  Double_t min=0,
1207  Double_t max=2.5);
1208  void DrawOne(Double_t cent=2.5,
1209  Double_t eta=0,
1210  Double_t ipz=0) const; //*MENU*
1215  ClassDef(AliTrackletDeltaWeights,1); // Base class of weights
1216 };
1217 
1218 //____________________________________________________________________
1220  const char* title)
1221  : AliTrackletBaseWeights(name, title),
1222  fHistos(),
1223  fCentAxis()
1224 {
1225  fHistos.SetOwner(true);
1226 }
1227 //____________________________________________________________________
1230  fHistos(),
1231  fCentAxis(o.fCentAxis)
1232 {
1233  fHistos.SetOwner(true);
1234  for (Int_t i = 0; i < o.fHistos.GetEntriesFast(); i++) {
1235  TH1* h = static_cast<TH1*>(o.fHistos.At(i));
1236  if (!h) continue;
1237  h = static_cast<TH1*>(h->Clone());
1238  h->SetDirectory(0);
1239  fHistos.Add(h);
1240  }
1241 }
1242 //____________________________________________________________________
1245 {
1246  if (&o == this) return *this;
1248  fCentAxis = o.fCentAxis;
1249  for (Int_t i = 0; i < o.fHistos.GetEntriesFast(); i++) {
1250  TH1* h = static_cast<TH1*>(o.fHistos.At(i));
1251  if (!h) continue;
1252  h = static_cast<TH1*>(h->Clone());
1253  h->SetDirectory(0);
1254  fHistos.Add(h);
1255  }
1256  return *this;
1257 }
1258 
1259 //____________________________________________________________________
1260 void
1262 {
1263  if (a.GetNbins() && a.GetXbins()->GetArray())
1264  SetCentAxis(a.GetNbins(), a.GetXbins()->GetArray());
1265  else
1266  SetCentAxis(a.GetNbins(), a.GetXmin(), a.GetXmax());
1267 }
1268 
1269 //____________________________________________________________________
1270 void
1272 {
1273  fCentAxis.Set(n, bins);
1274 }
1275 //____________________________________________________________________
1276 void
1278 {
1279  fCentAxis.Set(n, l, h);
1280 }
1281 
1282 //____________________________________________________________________
1283 Bool_t
1285 {
1286  if (bin < 1 || bin > fCentAxis.GetNbins()) {
1287  Warning("SetHisto", "Centrality bin %d out of range [%d,%d]",
1288  bin, 1, fCentAxis.GetNbins());
1289  return false;
1290  }
1291  Double_t c1 = fCentAxis.GetBinLowEdge(bin);
1292  Double_t c2 = fCentAxis.GetBinUpEdge(bin);
1293  TString name;
1294  name.Form("cent%03dd%02d_%03dd%02d",
1295  Int_t(c1), Int_t(c1 * 100) % 100,
1296  Int_t(c2), Int_t(c2 * 100) % 100);
1297  TH1* cpy = static_cast<TH1*>(h->Clone(name));
1298  cpy->SetDirectory(0);
1299  cpy->SetTitle(Form("%6.2f%% - %6.2f%%", c1, c2));
1300  fHistos.AddAtAndExpand(cpy, bin-1);
1301  return true;
1302 }
1303 
1304 //____________________________________________________________________
1305 TH1*
1307 {
1308  Int_t bin = FindBin(cent, &fCentAxis);
1309  if (bin >= fHistos.GetEntriesFast()) return 0;
1310  return static_cast<TH1*>(fHistos.At(bin-1));
1311 }
1312 
1313 //____________________________________________________________________
1314 Int_t
1316  const TAxis* axis)
1317 {
1318  if (!axis) return 1;
1319  Int_t bin = const_cast<TAxis*>(axis)->FindBin(value);
1320  if (bin < 1) bin = 1;
1321  else if (bin > axis->GetNbins()) bin = axis->GetNbins();
1322  return bin;
1323 }
1324 
1325 
1326 //____________________________________________________________________
1327 Double_t
1329  Double_t cent,
1330  Double_t ipZ,
1331  TH2* corr) const
1332 {
1333  if (tracklet->IsGenerated()) return 1;
1334  TH1* h = FindHisto(cent);
1335  if (!h) return 1;
1336 
1337  Double_t eta = tracklet->GetEta();
1338  Double_t delta = tracklet->GetDelta();
1339  Int_t etaBin = FindBin(eta, h->GetXaxis());
1340  Int_t deltaBin = FindBin(delta, h->GetYaxis());
1341  Int_t ipzBin = FindBin(ipZ, h->GetZaxis());
1342 
1343  return h->GetBinContent(etaBin, deltaBin, ipzBin);
1344 }
1345 
1346 //____________________________________________________________________
1347 TH1*
1348 AliTrackletDeltaWeights::Project(TH1* h, char which, const char* name)
1349 {
1350  TString n; n.Form("%s_%s",h->GetName(),name);
1351  TString t(h->GetTitle());
1352  TH1* ret = 0;
1353  TAxis* a0 = 0;
1354  TAxis* a1 = 0;
1355  TAxis* a2 = 0;
1356  typedef
1357  TH1D* (TH3::*ProjFunc)(const char*,Int_t,Int_t,Int_t,Int_t,Option_t*) const;
1358  ProjFunc func = 0;
1359  switch (which) {
1360  case 'x':
1361  a0 = h->GetXaxis();
1362  a1 = h->GetYaxis();
1363  a2 = h->GetZaxis();
1364  func = &TH3::ProjectionX;
1365  break;
1366  case 'y':
1367  a0 = h->GetYaxis();
1368  a1 = h->GetXaxis();
1369  a2 = h->GetZaxis();
1370  func = &TH3::ProjectionY;
1371  break;
1372  case 'z':
1373  a0 = h->GetZaxis();
1374  a1 = h->GetXaxis();
1375  a2 = h->GetYaxis();
1376  func = &TH3::ProjectionZ;
1377  break;
1378  }
1379  Int_t nx = a0->GetNbins();
1380  if (nx <= 1) return 0;
1381  if (a0->GetXbins() && a0->GetXbins()->GetArray())
1382  ret = new TH1D(n,t,nx,a0->GetXbins()->GetArray());
1383  else
1384  ret = new TH1D(n,t,nx,a0->GetXmin(),a0->GetXmax());
1385  ret->Sumw2();
1386  ret->SetXTitle(a0->GetTitle());
1387  static_cast<const TAttAxis*>(a0)->Copy(*(ret->GetXaxis()));
1388  ret->SetDirectory(0);
1389  ret->SetLineColor(h->GetMarkerColor());
1390  ret->SetMarkerColor(h->GetMarkerColor());
1391  ret->SetFillColor(kWhite);// color);
1392  ret->SetFillStyle(0);
1393  ret->SetMarkerStyle(20);
1394 
1395  for (Int_t i0 = 1; i0 <= ret->GetNbinsX(); i0++) {
1396  Int_t count = 0;
1397  Double_t sum = 0;
1398  Double_t sumE2 = 0;
1399  for (Int_t i1 = 1; i1 <= a1->GetNbins(); i1++) {
1400  for (Int_t i2 = 1; i2 <= a2->GetNbins(); i2++) {
1401  Int_t bin = 0;
1402  switch (which) {
1403  case 'x': bin = h->GetBin(i0, i1, i1); break;
1404  case 'y': bin = h->GetBin(i1, i0, i2); break;
1405  case 'z': bin = h->GetBin(i1, i2, i0); break;
1406  }
1407  Double_t c = h->GetBinContent(bin);
1408  Double_t e = h->GetBinError (bin);
1409  if (c < 1e-3 || e < 1e-6) continue;
1410  count++;
1411  sum += c;
1412  sumE2 += e*e;
1413  }
1414  }
1415  // ret->SetBinContent(i0, ret->GetBinContent(i0)/count);
1416  // ret->SetBinError (i0, ret->GetBinError (i0)/count);
1417  ret->SetBinContent(i0, sum/count);
1418  ret->SetBinError (i0, TMath::Sqrt(sumE2)/count);
1419  }
1420 
1421 
1422  // ret->Scale(1./(a1->GetNbins()*a2->GetNbins()));
1423  return ret;
1424 }
1425 
1426 //____________________________________________________________________
1427 void
1429 {
1431  gROOT->IncreaseDirLevel();
1432  gROOT->IndentLevel();
1433  Printf("%d centrality bins", fCentAxis.GetNbins());
1434  for (Int_t i = 0; i < fHistos.GetEntriesFast(); i++) {
1435  TH1* h = static_cast<TH1*>(fHistos.At(i));
1436  if (!h) continue;
1437  h->Print(option);
1438  }
1439  gROOT->DecreaseDirLevel();
1440 }
1441 //____________________________________________________________________
1442 void
1444  Double_t eta,
1445  Double_t ipz) const
1446 {
1447  TH1* h = FindHisto(cent);
1448  if (!h) return;
1449 
1450  Int_t etaBin = FindBin(eta, h->GetXaxis());
1451  Int_t ipzBin = FindBin(ipz, h->GetZaxis());
1452 
1453  TH1* tmp = 0;
1454  if (h->IsA()->InheritsFrom(TH2::Class()))
1455  tmp = static_cast<TH2*>(h)->ProjectionY("tmp", etaBin, etaBin);
1456  if (h->IsA()->InheritsFrom(TH3::Class()))
1457  tmp = static_cast<TH3*>(h)->ProjectionY("tmp", etaBin, etaBin);
1458  else
1459  tmp = static_cast<TH1*>(h->Clone("tmp"));
1460  if (!tmp) return;
1461  tmp->SetDirectory(0);
1462  tmp->Draw();
1463 }
1464 
1465 //____________________________________________________________________
1466 void
1468  Int_t nPad,
1469  THStack* stack,
1470  const char* xtitle,
1471  Double_t min,
1472  Double_t max)
1473 {
1474  if (!stack->GetHists()) return;
1475  stack->SetMinimum(min);
1476  stack->SetMaximum(max);
1477  pad->SetTicks();
1478  pad->SetGridy();
1479  pad->cd();
1480  stack->Draw("nostack");
1481  TH1* h = stack->GetHistogram();
1482  h->SetXTitle(xtitle);
1483  h->SetYTitle(stack->GetTitle());
1484  h->GetXaxis()->SetNdivisions(205);
1485  h->GetYaxis()->SetNdivisions(205);
1486  h->GetXaxis()->SetLabelOffset(0.005*pad->GetWNDC());
1487  h->GetXaxis()->SetTitleOffset(1.5*pad->GetWNDC());
1488  h->GetXaxis()->SetLabelSize(0.03/pad->GetWNDC()/*nPad*/);
1489  h->GetYaxis()->SetLabelSize(0.03/pad->GetWNDC()/*nPad*/);
1490  h->GetXaxis()->SetTitleSize(0.03/pad->GetWNDC()/*nPad*/);
1491  h->GetYaxis()->SetTitleSize(0.03/pad->GetWNDC()/*nPad*/);
1492  pad->Modified();
1493 }
1494 
1495 //____________________________________________________________________
1496 void
1498 {
1499  TVirtualPad* master = TVirtualPad::Pad();
1500  if (!master) {
1501  Warning("Draw", "No current pad to draw in");
1502  return;
1503  }
1504 
1505  TString opt(option); opt.ToUpper();
1506  THStack* stackEta = new THStack("eta", "#LTw#GT_{#Delta,IP_{z}}");
1507  THStack* stackDelta = new THStack("delta", "#LTw#GT_{#eta,IP_{z}}");
1508  THStack* stackIPz = new THStack("ipz", "#LTw#GT_{#eta,#Delta}");
1509 
1510  for (Int_t i = 0; i < fHistos.GetEntriesFast(); i++) {
1511  TH3* h = static_cast<TH3*>(fHistos.At(i));
1512  if (!h) continue;
1513 
1514  TH1* etaProj = Project(h, 'x', "eta");
1515  TH1* deltaProj = Project(h, 'y', "delta");
1516  TH1* ipzProj = Project(h, 'z', "ipz");
1517  stackEta ->Add(etaProj);
1518  stackDelta->Add(deltaProj);
1519  stackIPz ->Add(ipzProj);
1520  }
1521 
1522  Int_t nPad = 0;
1523  if (stackEta ->GetHists()) nPad++;
1524  if (stackDelta->GetHists()) nPad++;
1525  if (stackIPz ->GetHists()) nPad++;
1526 
1527  master->SetTopMargin(0.01);
1528  master->SetRightMargin(0.01);
1529  master->SetLeftMargin(0.07*nPad);
1530  if (opt.Contains("V"))
1531  master->Divide(1,nPad);
1532  else
1533  master->Divide(nPad,1, 0, 0);
1534 
1535  if (nPad >= 1) DrawStack(master->cd(1), nPad, stackEta, "#eta");
1536  if (nPad >= 2) DrawStack(master->cd(2), nPad, stackDelta, "#Delta");
1537  if (nPad >= 3) DrawStack(master->cd(3), nPad, stackIPz, "IP_{z}");
1538  master->GetPad(nPad)->SetRightMargin(0.01);
1539 
1540  master->Modified();
1541 }
1542 
1543 #endif
1544 //____________________________________________________________________
1545 //
1546 // EOF
1547 //
1548 
ClassDef(AliTrackletDeltaWeights, 1)
virtual Double_t CalcWeight(AliAODTracklet *tracklet, Double_t cent, Double_t ipz, TH2 *corr=0) const
Int_t pdg
void SetPtMode(UShort_t mode)
static Int_t FindBin(Double_t value, const TAxis *axis)
Real_t GetEta() const
double Double_t
Definition: External.C:58
Bool_t AddAbundanceWeight(Short_t pdg, const TH1D *h, UShort_t mode=0)
TCollection * Store(TCollection *out)
const char * title
Definition: MakeQAPdf.C:26
virtual TCollection * Store(TCollection *parent)
AliTrackletBaseWeights & operator=(const AliTrackletBaseWeights &o)
void Draw(Option_t *option="")
Definition: External.C:244
AliTrackletBaseWeights(const char *name, const char *title="Sim. tracklet weights")
void SetCalc(UChar_t mode=kProduct)
Bool_t AddStrangenessWeight(Short_t pdg, const TH1D *h, UShort_t mode=0)
void SetStrangenessMode(Short_t pdg, UShort_t mode)
void StoreMap(TCollection *parent, const char *name, PdgMap &m)
Bool_t SetHisto(Int_t bin, TH1 *h)
Bool_t RetrieveMap(TCollection *parent, const char *name, PdgMap &m)
void SetInverse(Bool_t inv)
TCanvas * c
Definition: TestFitELoss.C:172
virtual Double_t CalcWeight(AliAODTracklet *tracklet, Double_t cent, Double_t ipz, TH2 *corr) const =0
void DrawStack(TVirtualPad *pad, Int_t nPad, THStack *stack, const char *xtitle, Double_t min=0, Double_t max=2.5)
AliTrackletDeltaWeights & operator=(const AliTrackletDeltaWeights &o)
ClassDef(AliTrackletPtPidStrWeights, 3)
void Print(Option_t *option="") const
void Draw(Option_t *option="")
virtual Real_t GetParentPt(Bool_t second=false) const
void SetMask(UChar_t mask)
UShort_t T(UShort_t m, UShort_t t)
Definition: RingBits.C:60
virtual TCollection * Retrieve(TCollection *in)
Bool_t IsSimulated() const
void ModStack(THStack *stack)
Tracklet AOD object.
void SetCentAxis(const TAxis &axis)
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
TCollection * Retrieve(TCollection *in)
AliTrackletBaseWeights(const AliTrackletBaseWeights &o)
Definition: External.C:228
Definition: External.C:212
Bool_t IsMeasured() const
virtual Short_t GetParentPdg(Bool_t second=false) const
virtual void Print(Option_t *option="") const
virtual Double_t CalcWeight(AliAODTracklet *tracklet, Double_t cent, Double_t ipZ, TH2 *corr=0) const
Int_t mode
Definition: anaM.C:40
Bool_t SetPtWeight(const TH2D *h, UShort_t mode=0)
Double_t GetStrangenessWeight(UShort_t apdg, Double_t cent) const
TH1D * GetPdgHist(const PdgMap &m, Short_t pdg) const
void SetVeto(UChar_t veto)
Bool_t AddPdgWeight(PdgMap &map, Short_t pdg, const TH1D *w, UShort_t mode=0)
short Short_t
Definition: External.C:23
Int_t colors[nPtBins]
Double_t LookupWeight(AliAODTracklet *tracklet, Double_t cent, Double_t ipz, TH2 *corr=0) const
UChar_t GetFlags() const
Bool_t IsGenerated() const
void SetPdgMode(PdgMap &map, Short_t pdg, UShort_t mode)
std::map< short, TH1D * > PdgMap
void DrawOne(Double_t cent=2.5, Double_t eta=0, Double_t ipz=0) const
Definition: External.C:220
void PrintHist(const char *tag, TH1 *h) const
void SetAbundanceMode(Short_t pdg, UShort_t mode)
void Print(Option_t *option="") const
virtual Double_t GetPdgWeight(const PdgMap &m, UShort_t apdg, Double_t cent) const
unsigned short UShort_t
Definition: External.C:28
const char Option_t
Definition: External.C:48
void PrintMap(const PdgMap &m, const char *name, Option_t *options="") const
Real_t GetDelta() const
TList * GetHists(UShort_t flags, const char *var, const char *stackName="result", const char *sub="")
Definition: Compare.C:81
bool Bool_t
Definition: External.C:53
TH1 * FindHisto(Double_t cent) const
Double_t GetAbundanceWeight(UShort_t apdg, Double_t cent) const
AliTrackletPtPidStrWeights & operator=(const AliTrackletPtPidStrWeights &o)
Bool_t CheckTracklet(const AliAODTracklet *tracklet) const
TH1 * Project(TH1 *h, char which, const char *name)
Definition: External.C:196
ClassDef(AliTrackletBaseWeights, 1)