AliPhysics  026afea (026afea)
 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 <TList.h>
19 # include <TBrowser.h>
20 # include <TCanvas.h>
21 # include <THStack.h>
22 # include <TROOT.h>
23 # include <TClass.h>
24 # include "AliAODTracklet.C"
25 # include <TParameter.h>
26 #else
27 class TH1D;
28 class TH2D;
29 class TH2;
30 class TList;
31 class TBrowser;
32 class THStack;
33 class AliAODTracklet;
34 class TCanvas; // Autoload
35 #endif
36 
42 class AliTrackletWeights : public TNamed
43 {
44 public:
48  enum EMode {
50  kUp = (0x1) << 14,
52  kDown = (0x2) << 14,
54  kDisabled = (0x4) << 14
55  };
59  enum ECalc {
102  };
103 
105  typedef std::map<short,TH1D*> PdgMap;
111  : TNamed(),
112  fPt(0),
113  fAbundance(),
114  fStrangeness(),
115  fCalc(kProduct)
116  {}
123  AliTrackletWeights(const char* name,
124  const char* title="Sim. tracklet weights");
134  virtual ~AliTrackletWeights() {}
143 
153  virtual Double_t LookupWeight(AliAODTracklet* tracklet,
154  Double_t cent,
155  TH2* corr=0) const;
166  {
167  AddPdgWeight(fAbundance, pdg, h, mode);
168  }
179  {
180  AddPdgWeight(fStrangeness, pdg, h, mode);
181  }
190  Bool_t SetPtWeight(const TH2D* h, UShort_t mode=0);
204  {
205  SetPdgMode(fAbundance, pdg, mode);
206  }
219  {
220  SetPdgMode(fStrangeness, pdg, mode);
221  }
232  void SetPtMode(UShort_t mode);
280  void SetCalc(UChar_t mode=kProduct) { fCalc = mode; }
286  void SetMask(UChar_t mask) { fMask = mask; }
292  void SetVeto(UChar_t veto) { fVeto = veto; }
298  void Draw(Option_t* option=""); //*MENU*
304  void Print(Option_t* option="") const; //*MENU*
308  Bool_t IsFolder() const { return true; }
314  void Store(TCollection* out);
333  {
334  return GetPdgWeight(fAbundance, apdg, cent);
335  }
346  {
347  return GetPdgWeight(fStrangeness, apdg, cent);
348  }
349 
350 protected:
351 
362  virtual Double_t LookupWeight(Double_t pT, Short_t pdg, Double_t cent) const;
371  TH1D* GetPdgHist(const PdgMap& m, Short_t pdg) const;
381  virtual Double_t GetPdgWeight(const PdgMap& m,
382  UShort_t apdg,
383  Double_t cent) const;
394  Bool_t AddPdgWeight(PdgMap& map, Short_t pdg, const TH1D* w, UShort_t mode=0);
406  void SetPdgMode(PdgMap& map, Short_t pdg, UShort_t mode);
414  void StoreMap(TCollection* parent, const char* name, PdgMap& m);
424  Bool_t RetrieveMap(TCollection* parent, const char* name, PdgMap& m);
430  void ModStack(THStack* stack);
438  void PrintMap(const PdgMap& m, const char* name, Option_t* options="") const;
445  void PrintHist(const char* tag, TH1* h) const;
446 
447  TH2D* fPt; // Weight by centrality and pT
448  PdgMap fAbundance; // Map for abundance weight
449  PdgMap fStrangeness; // Map for strangeness weight
450  UChar_t fCalc; // Whether the square of the weight is calculated
451  UChar_t fMask; // Which particles to take
452  UChar_t fVeto; // Which particles not to take
453  ClassDef(AliTrackletWeights,3); // Weighs for tracklet analysis
454 };
455 
456 //____________________________________________________________________
458  const char* title)
459  : TNamed(name, title),
460  fPt(0),
461  fAbundance(),
462  fStrangeness(),
463  fCalc(kProduct),
464  fMask(0xFF),
465  fVeto(0x0)
466 {}
467 
468 //____________________________________________________________________
470  : TNamed(o),
471  fPt(0),
472  fAbundance(),
473  fStrangeness(),
474  fCalc(o.fCalc),
475  fMask(o.fMask),
476  fVeto(o.fVeto)
477 {
478  UInt_t mask = kUp|kDown|kDisabled;
479  SetPtWeight(o.fPt, o.fPt->TestBits(mask));
480  for (PdgMap::const_iterator i = o.fAbundance.begin();
481  i != o.fAbundance.end(); ++i)
482  AddAbundanceWeight(i->first, i->second, i->second->TestBits(mask));
483  for (PdgMap::const_iterator i = o.fStrangeness.begin();
484  i != o.fStrangeness.end(); ++i)
485  AddStrangenessWeight(i->first, i->second, i->second->TestBits(mask));
486 }
487 //____________________________________________________________________
489 {
490  if (&o == this) return *this;
491  fName = o.fName;
492  fTitle = o.fTitle;
493  fCalc = o.fCalc;
494  fMask = o.fMask;
495  fVeto = o.fVeto;
496  UInt_t mask = kUp|kDown|kDisabled;
497  SetPtWeight(o.fPt, o.fPt->TestBits(mask));
498  for (PdgMap::const_iterator i = o.fAbundance.begin();
499  i != o.fAbundance.end(); ++i)
500  AddAbundanceWeight(i->first, i->second, i->second->TestBits(mask));
501  for (PdgMap::const_iterator i = o.fStrangeness.begin();
502  i != o.fStrangeness.end(); ++i)
503  AddStrangenessWeight(i->first, i->second, i->second->TestBits(mask));
504  return *this;
505 }
506 //____________________________________________________________________
508  Short_t pdg,
509  const TH1D* w,
510  UShort_t mode)
511 {
512  UShort_t apdg = TMath::Abs(pdg);
513  TH1D* copy = static_cast<TH1D*>(w->Clone(Form("w%d", apdg)));
514  copy->SetDirectory(0);
515  copy->SetXTitle("Centrality [%]");
516  copy->SetBit(mode);
517  PdgMap::iterator i = m.find(apdg);
518  if (i != m.end()) {
519  Warning("AddPdgWeight", "Replacing weight for %d", pdg);
520  delete i->second;
521  i->second = 0;
522  m.erase(i);
523  }
524  m[apdg] = copy;
525  return true;
526 }
527 
528 //____________________________________________________________________
530 {
531  if (fPt) {
532  delete fPt;
533  fPt = 0;
534  }
535  if (!h) return false;
536  fPt = static_cast<TH2D*>(h->Clone("centPt"));
537  fPt->SetDirectory(0);
538  fPt->SetXTitle("Centrality [%]");
539  fPt->SetYTitle("#it{p}_{T} [GeV]");
540  fPt->SetBit(mode);
541  return true;
542 }
543 //____________________________________________________________________
545 {
546  if (fPt) fPt->SetBit(mode);
547 }
548 //____________________________________________________________________
550  Short_t pdg,
551  UShort_t mode)
552 {
553  TH1* h = GetPdgHist(m, pdg);
554  if (h) h->SetBit(mode);
555 }
556 //____________________________________________________________________
558 {
559  UShort_t apdg = TMath::Abs(pdg);
560  PdgMap::const_iterator i = m.find(apdg);
561  if (i == m.end()) return 0;
562  return i->second;
563 }
564 
565 //____________________________________________________________________
567  UShort_t apdg,
568  Double_t cent) const
569 {
570  if (m.size() < 1) return 1;
571  TH1D* h = GetPdgHist(m, apdg);
572  if (!h || h->TestBit(kDisabled)) return 1;
573 
574 #if 1
575  Int_t bC = h->GetXaxis()->FindBin(cent);
576  if (bC < 1 || bC > h->GetNbinsX()) return 1;
577 #else
578  Int_t bC = 1;
579 #endif
580  Double_t add = (h->TestBit(kUp) ? +1 :
581  h->TestBit(kDown) ? -1 : 0) * h->GetBinError(bC);
582  return h->GetBinContent(bC) + add;
583 }
584 //____________________________________________________________________
585 Double_t
587 {
588  Double_t w = 1;
589  if (fPt && !fPt->TestBit(kDisabled)) {
590  Int_t bC = fPt->GetXaxis()->FindBin(cent);
591  Int_t bpT = fPt->GetYaxis()->FindBin(pT);
592  if (bC >= 1 && bC <= fPt->GetXaxis()->GetNbins() &&
593  bpT >= 1 && bpT <= fPt->GetYaxis()->GetNbins()) {
594  Double_t fac = (pT >= 0.05 ? 1 :
595  (fPt->TestBit(kUp) ? 1.3 :
596  fPt->TestBit(kDown) ? 0.7 : 1));
597  w *= fac*fPt->GetBinContent(bC, bpT);
598  }
599  }
600  UShort_t apdg = TMath::Abs(pdg);
601 
602  w *= GetPdgWeight(fAbundance, apdg, cent);
603  w *= GetPdgWeight(fStrangeness, apdg, cent);
604  // Printf("Weight of pT=%6.3f pdg=%5d cent=%5.1f -> %f", pT, pdg, cent, w);
605  return w;
606 }
607 
608 //____________________________________________________________________
609 Double_t
611  Double_t cent,
612  TH2* corr) const
613 {
614 #if 0
615  if (!tracklet->IsSimulated()) {
616  Warning("LookupWeight", "Not a simulated tracklet");
617  return 1;
618  }
619 #endif
620  UChar_t flags = tracklet->GetFlags();
621  if (fMask != 0xFF && (fMask & flags) == 0) {
622  // Info("LookupWeight", "Tracklet 0x%02x does not fullfill mask 0x%02x",
623  // flags, fMask);
624  return 1;
625  }
626  if ((fVeto & flags) != 0) {
627  // Info("LookupWeight", "Tracklet 0x%02x vetoed by 0x%02x",flags, fVeto);
628  return 1;
629  }
630 
631  Double_t w1 = 1, w2 = 1;
632 
633  // AliAODMCTracklet* mc = static_cast<AliAODMCTracklet*>(tracklet);
634  AliAODTracklet* mc = tracklet;
635  Short_t pdg1 = mc->GetParentPdg();
636  Short_t pdg2 = mc->GetParentPdg(true);
637  if (pdg1>0) w1 = LookupWeight(mc->GetParentPt(), pdg1,cent);
638  if (pdg2>0) w2 = LookupWeight(mc->GetParentPt(true),pdg2,cent);
639  else if (fCalc!=kProduct) w2 = w1;
640 
641  if (corr && mc->IsMeasured()) corr->Fill(w1, w2);
642 
643  switch (fCalc) {
644  case kProduct: return w1 * w2;
645  case kSquare: return TMath::Sqrt(w1 * w2);
646  case kSum: return 1+(w1-1)+(w2-1);
647  case kAverage: return 1+((w1-1)+(w2-1))/2;
648  }
649  return 1;
650 }
651 
652 //____________________________________________________________________
653 void AliTrackletWeights::ModStack(THStack* stack)
654 {
655  if (!stack || !stack->GetHists()) return;
656  TIter next(stack->GetHists());
657  TH1* hist = 0;
658  Color_t colors[] = { kPink+2, kRed+2, kOrange+2, kYellow+2,
659  kSpring+2, kGreen+2, kTeal+2, kCyan+2,
660  kBlue+2, kViolet+2 };
661  Int_t idx = 0;
662  while ((hist = static_cast<TH1*>(next()))) {
663  Color_t c = colors[idx % 10];
664  idx++;
665  hist->SetLineColor(c);
666  hist->SetMarkerColor(c);
667  hist->SetMarkerStyle(20+idx%20);
668  }
669 }
670 
671 //____________________________________________________________________
672 void
674 {
675  TList* top = new TList;
676  top->SetName(GetName());
677  top->SetOwner(true);
678  parent->Add(top);
679  top->Add(new TParameter<int>("calc", fCalc, 'f'));
680  top->Add(new TParameter<int>("mask", fMask, 'f'));
681  top->Add(new TParameter<int>("veto", fVeto, 'f'));
682 
683  if (fPt) {
684  TH2* copy = static_cast<TH2*>(fPt->Clone("centPt"));
685  copy->SetDirectory(0);
686  copy->SetBinContent(0,0,1); // For counting merges
687  top->Add(copy);
688  }
689  StoreMap(top, "abundance", fAbundance);
690  StoreMap(top, "strangeness", fStrangeness);
691 }
692 
693 //____________________________________________________________________
694 void
695 AliTrackletWeights::StoreMap(TCollection* parent, const char* name, PdgMap& m)
696 {
697  TList* top = new TList;
698  top->SetName(name);
699  top->SetOwner(true);
700  parent->Add(top);
701  for (PdgMap::const_iterator i = m.begin(); i != m.end(); ++i) {
702  TH1* copy = static_cast<TH1*>(i->second->Clone(Form("w%d",i->first)));
703  copy->SetDirectory(0);
704  copy->SetBinContent(0,1); // For counting merges
705  top->Add(copy);
706  }
707 }
708 
709 #ifndef __CINT__
710 namespace {
711  template <typename T>
712  T GetPar(const char* name, TList* l)
713  {
714  TObject* o = l->FindObject(name);
715  if (!o) {
716  Warning("GetPar", "Didn't find parameter %s in %s",
717  name, l->GetName());
718  return T();
719  }
720  TClass* cls = TParameter<T>::Class();
721  if (!o->IsA()->InheritsFrom(cls)) {
722  Warning("GetPar", "Object %s is a %s, not a %s",
723  name, o->ClassName(), cls->GetName());
724  return T();
725  }
726  TParameter<T>* p = static_cast<TParameter<T>*>(o);
727  return p->GetVal();
728  }
729 }
730 #endif
731 
732 //____________________________________________________________________
733 Bool_t
735 {
736  TList* top = static_cast<TList*>(parent->FindObject(GetName()));
737  if (!top) {
738  Warning("Retrieve",
739  "Collection %s not found in %s", GetName(), parent->GetName());
740  parent->ls();
741  return false;
742  }
743  fCalc = GetPar<int>("calc", top);
744  fMask = GetPar<int>("mask", top);
745  fVeto = GetPar<int>("veto", top);
746 
747  fPt = static_cast<TH2D*>(top->FindObject("centPt"));
748  if (!fPt)
749  Warning("Retrieve","centPt histogram not found in %s", GetName());
750  else {
751  fPt->SetDirectory(0);
752  Double_t scale = fPt->GetBinContent(0,0);
753  fPt->Scale(1/scale); // Counting merges
754  // fPt->SetBinContent(0,0,1); // Zero merger count
755  }
756  if (!RetrieveMap(top, "abundance", fAbundance)) return false;
757  if (!RetrieveMap(top, "strangeness", fStrangeness)) return false;
758 
759  return true;
760 }
761 //____________________________________________________________________
762 Bool_t
764  const char* name,
765  PdgMap& m)
766 {
767  m.clear();
768  TList* top = static_cast<TList*>(parent->FindObject(name));
769  if (!top) {
770  Warning("RetrieveMap",
771  "Collection %s not found in %s", name, parent->GetName());
772  return true;
773  }
774  TIter next(top);
775  TObject* o = 0;
776  while ((o = next())) {
777  if (!o->IsA()->InheritsFrom(TH1D::Class())) continue;
778  TH1D* copy = static_cast<TH1D*>(o);
779  copy->SetDirectory(0);
780  Double_t scale = copy->GetBinContent(0);
781  copy->Scale(1/scale); // Counting merges
782  // copy->SetBinContent(0,1); // Zero merger count
783  TString nme(copy->GetName());
784  nme.Remove(0,1);
785  Int_t apdg = nme.Atoi();
786  AddPdgWeight(m, apdg, copy, copy->TestBits(kUp|kDown|kDisabled));
787  }
788  return true;
789 }
790 
791 //____________________________________________________________________
792 void
794 {
795  TVirtualPad* master = TVirtualPad::Pad();
796  if (!master) {
797  Warning("Draw", "No current pad to draw in");
798  return;
799  }
800 
801  Int_t nPad = 3;
802  Int_t iPad = 0;
803  if (fAbundance.size() <= 0) nPad--;
804  if (fStrangeness.size() <= 0) nPad--;
805  TString opt(option); opt.ToUpper();
806 
807  master->Divide(nPad,1);
808  if (fPt) {
809  THStack* stack = new THStack(fPt, (opt.Contains("C") ? "x" : "y"));
810  ModStack(stack);
811  stack->SetTitle("#it{p}_{T} weights");
812  master->cd(++iPad);
813  stack->Draw("nostack");
814  stack->GetHistogram()->SetXTitle((opt.Contains("C") ?
815  "Centrality [%]" : "#it{p}_{T}"));
816  master->GetPad(iPad)->BuildLegend();
817  master->GetPad(iPad)->Modified();
818  }
819 
820  if (fAbundance.size() > 0) {
821  THStack* ha = new THStack("abundance", "Abundance weights");
822  for (PdgMap::const_iterator i = fAbundance.begin();
823  i != fAbundance.end(); ++i) {
824  ha->Add(i->second);
825  }
826  ModStack(ha);
827  master->cd(++iPad);
828  ha->Draw("nostack");
829  ha->GetHistogram()->SetXTitle("Centrality [%]");
830  master->GetPad(iPad)->BuildLegend();
831  }
832 
833  if (fStrangeness.size() > 0) {
834  THStack* hs = new THStack("strangeness", "Strangeness weights");
835  for (PdgMap::const_iterator i = fStrangeness.begin();
836  i != fStrangeness.end(); ++i) {
837  hs->Add(i->second);
838  }
839  ModStack(hs);
840  master->cd(++iPad);
841  hs->Draw("nostack");
842  hs->GetHistogram()->SetXTitle("Centrality [%]");
843  master->GetPad(iPad)->BuildLegend();
844  }
845  master->Modified();
846 }
847 
848 //____________________________________________________________________
849 void
851 {
852  gROOT->IndentLevel();
853  Printf("%s : %s", ClassName(), GetName());
854  Printf(" Weight calculation: %s",
855  fCalc == kProduct ? "product" :
856  fCalc == kSquare ? "square" :
857  fCalc == kSum ? "sum" : "average");
858  Printf(" Tracklet mask: 0x%02x", fMask);
859  Printf(" Tracklet veto: 0x%02x", fVeto);
860  gROOT->IncreaseDirLevel();
861  PrintHist("pT", fPt);
862  PrintMap(fAbundance, "Abundance", option);
863  PrintMap(fStrangeness, "Strangeness", option);
864 
865  gROOT->DecreaseDirLevel();
866 }
867 //____________________________________________________________________
868 void
869 AliTrackletWeights::PrintHist(const char* tag, TH1* h) const
870 {
871  if (!h) return;
872  gROOT->IndentLevel();
873  Printf("%10s (%c): %p %s/%s", tag,
874  h->TestBit(kDisabled) ? '0' :
875  h->TestBit(kUp) ? '+' :
876  h->TestBit(kDown) ? '-' : '=',
877  h, h->GetName(), h->GetTitle());
878 
879 }
880 //____________________________________________________________________
881 void
883  const char* name,
884  Option_t* option) const
885 {
886  gROOT->IndentLevel();
887  Printf("Map of PDG codes: %s", name);
888  gROOT->IncreaseDirLevel();
889  for (PdgMap::const_iterator i = m.begin(); i != m.end(); ++i)
890  PrintHist(Form("%10d", i->first), i->second);
891 
892  gROOT->DecreaseDirLevel();
893 }
894 #endif
895 //____________________________________________________________________
896 //
897 // EOF
898 //
899 
Bool_t SetPtWeight(const TH2D *h, UShort_t mode=0)
Bool_t IsFolder() const
Int_t pdg
double Double_t
Definition: External.C:58
Double_t GetAbundanceWeight(UShort_t apdg, Double_t cent) const
void SetMask(UChar_t mask)
const char * title
Definition: MakeQAPdf.C:26
Bool_t RetrieveMap(TCollection *parent, const char *name, PdgMap &m)
void Draw(Option_t *option="")
TCanvas * c
Definition: TestFitELoss.C:172
void SetPtMode(UShort_t mode)
void PrintMap(const PdgMap &m, const char *name, Option_t *options="") const
void StoreMap(TCollection *parent, const char *name, PdgMap &m)
virtual Real_t GetParentPt(Bool_t second=false) const
UShort_t T(UShort_t m, UShort_t t)
Definition: RingBits.C:60
TH1D * GetPdgHist(const PdgMap &m, Short_t pdg) const
Bool_t IsSimulated() const
Tracklet AOD object.
int Int_t
Definition: External.C:63
void SetPdgMode(PdgMap &map, Short_t pdg, UShort_t mode)
virtual Double_t GetPdgWeight(const PdgMap &m, UShort_t apdg, Double_t cent) const
unsigned int UInt_t
Definition: External.C:33
std::map< short, TH1D * > PdgMap
Double_t GetStrangenessWeight(UShort_t apdg, Double_t cent) const
void ModStack(THStack *stack)
Definition: External.C:228
Definition: External.C:212
void SetCalc(UChar_t mode=kProduct)
Bool_t Retrieve(TCollection *in)
Bool_t IsMeasured() const
virtual Short_t GetParentPdg(Bool_t second=false) const
Int_t mode
Definition: anaM.C:40
virtual Double_t LookupWeight(AliAODTracklet *tracklet, Double_t cent, TH2 *corr=0) const
short Short_t
Definition: External.C:23
UChar_t GetFlags() const
Bool_t AddPdgWeight(PdgMap &map, Short_t pdg, const TH1D *w, UShort_t mode=0)
void Store(TCollection *out)
void SetStrangenessMode(Short_t pdg, UShort_t mode)
Bool_t AddStrangenessWeight(Short_t pdg, const TH1D *h, UShort_t mode=0)
Definition: External.C:220
void SetAbundanceMode(Short_t pdg, UShort_t mode)
void SetVeto(UChar_t veto)
unsigned short UShort_t
Definition: External.C:28
void PrintHist(const char *tag, TH1 *h) const
const char Option_t
Definition: External.C:48
Bool_t AddAbundanceWeight(Short_t pdg, const TH1D *h, UShort_t mode=0)
bool Bool_t
Definition: External.C:53
void Print(Option_t *option="") const
ClassDef(AliTrackletWeights, 3)
Definition: External.C:196
AliTrackletWeights & operator=(const AliTrackletWeights &o)