AliPhysics  781d0c7 (781d0c7)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliTrackletAODdNdeta.C
Go to the documentation of this file.
1 
12 #include <AliAnalysisTaskSE.h>
13 #include <AliTrackletAODUtils.C>
14 #include <AliAODTracklet.C>
15 #ifndef __CINT__
16 #include <cctype>
17 #include "AliAODSimpleHeader.C"
18 #include <AliVVertex.h>
19 #include <AliVertex.h>
20 #include <AliAnalysisManager.h>
21 #include <AliVEventHandler.h>
22 #include <AliInputEventHandler.h>
23 #include <AliMultSelection.h>
24 #include <AliCentrality.h>
25 #include <AliLog.h>
26 #include <TClonesArray.h>
27 #include "AliTrackletWeights.C"
28 #include <TUrl.h>
29 #include <TFile.h>
30 #include <TGraphErrors.h>
31 #else
32 // class AliAODTracklet;
34 class AliVEvent;
35 class AliMultSelection; // Auto-load
36 class TClonesArray;
37 #endif
38 
39 //====================================================================
46  public AliTrackletAODUtils
47 {
48 public:
52  enum {
62  kPrimaryMask = 0x0,
71  };
72 
74  typedef TList Container;
82  AliTrackletAODdNdeta(const char* name);
92  virtual ~AliTrackletAODdNdeta() {}
110  virtual void Print(Option_t* option="") const;
117  Bool_t Connect(const char* sumFile=0,const char* resFile=0);
128  static AliTrackletAODdNdeta* Create(Bool_t mc=false,
129  const char* weights=0,
130  const char* sumFile=0,
131  const char* resFile=0);
132  /* @} */
133 
134  // -----------------------------------------------------------------
147  void UserExec(Option_t*);
152  void FinishTaskOutput() { /*WorkerFinalize();*/ }
157  void Terminate(Option_t*);
158  /* @} */
159 
160  // -----------------------------------------------------------------
170  void SetCentralityMethod(const TString& name) { fCentMethod = name; }
178  {
179  SetAxis(fCentAxis,n,bins);
180  }
186  void SetCentralityAxis(const TString& spec)
187  {
188  SetAxis(fCentAxis, spec, "-:,");
189  }
197  void SetIPzAxis(Int_t n, Double_t min, Double_t max)
198  {
199  SetAxis(fIPzAxis, n, min, max);
200  }
207  void SetIPzAxis(Int_t n, Double_t max)
208  {
209  SetAxis(fIPzAxis, n, max);
210  }
216  void SetIPzAxis(const TString& spec)
217  {
218  SetAxis(fIPzAxis, spec);
219  }
227  void SetPhiAxis(Int_t n, Double_t min, Double_t max)
228  {
229  SetAxis(fPhiAxis, n, min, max);
230  }
237  void SetPhiAxis(Int_t n, Double_t max)
238  {
239  SetAxis(fPhiAxis, n, max);
240  }
248  void SetEtaAxis(Int_t n, Double_t min, Double_t max)
249  {
250  SetAxis(fEtaAxis, n, min, max);
251  }
258  void SetEtaAxis(Int_t n, Double_t max)
259  {
260  SetAxis(fEtaAxis, n, max);
261  }
267  void SetEtaAxis(const TString& spec)
268  {
269  SetAxis(fEtaAxis, spec);
270  }
276  void SetDPhiShift(Double_t x=0.0045) { fDPhiShift = x; }
288  void SetDeltaCut(Double_t x=1.5) { fDeltaCut = x; }
294  void SetMaxDelta(Double_t x=25) { fMaxDelta = x; }
300  void SetTailDelta(Double_t x=5) { fTailDelta = x; }
306  void SetTailMaximum(Double_t x=-1) { fTailMax = x; }
320  void SetAbsMinCent(Double_t x=-1) { fAbsMinCent = x; }
326  virtual void SetWeights(AliTrackletBaseWeights* w) {};
332  virtual void SetWeightCalc(UChar_t mode=0) {}
338  virtual void SetWeightMask(UChar_t mask=0xFF) {}
346  virtual void SetWeightInverse(Bool_t inv) {}
352  virtual void SetWeightVeto(UChar_t veto=0xFF) {}
353  /* @} */
354  //__________________________________________________________________
358  struct Sub : public TObject
359  {
364  Sub(const char* name="") : fName(name), fContainer(0), fDebug(0) {}
370  Sub(const Sub& o) : TObject(o), fName(o.fName), fContainer(0), fDebug(0) {}
374  virtual ~Sub() {}
380  Sub& operator=(const Sub&) { return *this; }
384  const char* GetName() const { return fName.Data(); }
395  virtual Bool_t WorkerInit(Container* parent,
396  const TAxis& etaAxis,
397  const TAxis& ipzAxis,
398  const TAxis& deltaAxis)
399  {
400  fContainer = new Container;
401  fContainer->SetName(fName);
402  fContainer->SetOwner();
403  if (parent) parent->Add(fContainer);
404  return true;
405  }
416  virtual Bool_t ProcessTracklet(AliAODTracklet* tracklet,
417  Double_t ipz,
418  UShort_t signal,
419  Double_t weight) = 0;
429  virtual Bool_t FinalizeInit(Container* parent) = 0;
440  virtual Bool_t MasterFinalize(Container* parent,
441  TH1* ipz,
442  Double_t tailCut,
443  Double_t tailMax) = 0;
449  virtual void SetDebug(UShort_t lvl) { fDebug = lvl; }
450  protected:
457  ClassDef(Sub,1);
458  };
459 
460  //__________________________________________________________________
470  struct Histos : public Sub
471  {
476  Histos(const char* name="", Color_t color=kBlack, Style_t style=1,
477  UChar_t mask=0, UChar_t veto=0)
478  : Sub(name),
479  fMask(mask),
480  fVeto(veto),
481  fEtaIPz(0),
482  fEtaDeltaIPz(0),
483  fEtaDeltaPdg(0),
484  fEtaPdgIPz(0),
485  fEtaPdg(0),
486  fEtaPt(0),
487  fColor(color),
488  fStyle(style)
489  {}
495  Histos(const Histos& o)
496  : Sub(o),
497  fMask(o.fMask),
498  fVeto(o.fVeto),
499  fEtaIPz(0),
500  fEtaDeltaIPz(0),
501  fEtaDeltaPdg(0),
502  fEtaPdgIPz(0),
503  fEtaPdg(0),
504  fEtaPt(0),
505  fColor(o.fColor),
506  fStyle(o.fStyle)
507  {}
511  virtual ~Histos() {}
517  Histos& operator=(const Histos&) { return *this; }
522  {
523  return fMask == kMeasuredMask && fVeto == kMeasuredVeto;
524  }
529  {
530  return fMask == kInjectedMask && fVeto == kInjectedVeto;
531  }
536  {
538  }
544  {
545  return fMask == kDistinctMask && fVeto == kDistinctVeto;
546  }
551  {
552  return fMask == kPrimaryMask && fVeto == kPrimaryVeto;
553  }
558  {
559  return fMask == kSecondaryMask && fVeto == kSecondaryVeto;
560  }
565  {
566  return fMask == kGeneratedMask && fVeto == kGeneratedVeto;
567  }
578  Bool_t WorkerInit(Container* parent,
579  const TAxis& etaAxis,
580  const TAxis& ipzAxis,
581  const TAxis& deltaAxis);
588  void SetAttr(Color_t c, Style_t m);
600  Double_t ipz,
601  UShort_t signal,
602  Double_t weight);
612  Bool_t FinalizeInit(Container* parent);
624  TH1* ipz,
625  Double_t tailCut,
626  Double_t tailMax);
627 
628  UChar_t GetMask() const { return fMask; }
629  UChar_t GetVeto() const { return fVeto; }
635  void Print(Option_t* option="") const;
636  protected:
639  Int_t nEvents,
640  Double_t tailDelta,
641  Double_t tailMax,
642  const TString& pre,
643  const TString& tit);
645  Int_t nEvents,
646  Double_t tailCut,
647  Double_t tailMax);
649  TH1* ipz,
650  const TString& shn);
651  const UChar_t fMask;
652  const UChar_t fVeto;
659  public:
660  Color_t fColor;
661  Style_t fStyle;
662 
663  ClassDef(Histos,1);
664  };
665 
666 
667  //__________________________________________________________________
683  struct CentBin : public Sub
684  {
689  : Sub(""),
690  fSubs(0),
691  fLow(0),
692  fHigh(0),
693  fIPz(0),
694  fCent(0),
695  fCentIPz(0),
696  fMeasured(0),
697  fInjection(0)
698  {
699  }
706  CentBin(Double_t c1, Double_t c2);
712  CentBin(const CentBin& o)
713  : Sub(o),
714  fSubs(0),
715  fLow(o.fLow),
716  fIPz(0),
717  fCent(0),
718  fCentIPz(0),
719  fHigh(o.fHigh),
720  fMeasured(0),
721  fInjection(0)
722  {}
726  virtual ~CentBin() {}
732  CentBin& operator=(const CentBin&) { return *this; }
744  Histos* MakeHistos(const char* name,
745  Color_t color,
746  Style_t style,
747  UShort_t mask,
748  UShort_t veto);
759  Bool_t WorkerInit(Container* parent,
760  const TAxis& etaAxis,
761  const TAxis& ipzAxis,
762  const TAxis& deltaAxis);
768  Bool_t IsAllBin() const;
777  Bool_t Accept(Double_t cent, Double_t ipz);
789  Double_t ipz,
790  UShort_t signal,
791  Double_t weight);
801  Bool_t FinalizeInit(Container* parent);
813  TH1* ipz,
814  Double_t tailCut,
815  Double_t tailMax);
829  Container* measCont,
830  Container* genCont,
831  Histos* h,
832  Double_t tailCut,
833  Double_t tailMax);
839  void Print(Option_t* option="") const;
845  virtual void SetDebug(UShort_t lvl);
846  protected:
852  TProfile* fCentIPz;
855 
856  ClassDef(CentBin,2);
857  };
858 
859 protected:
860  // -----------------------------------------------------------------
870  virtual Bool_t WorkerInit();
877  Bool_t InitCentBins(Container* existing);
887  {
888  return new CentBin(c1, c2);
889  }
901  virtual CentBin* InitCentBin(Int_t bin,
902  Float_t c1,
903  Float_t c2,
904  Container* existing,
905  const TAxis& deltaAxis);
906  /* @} */
907 
908  // -----------------------------------------------------------------
913  virtual const char* GetBranchName() const { return "AliAODTracklets"; }
923  Bool_t CheckEvent(Double_t& cent,
924  const AliVVertex*& ip,
925  TClonesArray*& tracklets);
933  TClonesArray* FindTracklets(AliVEvent* event);
943  const AliVVertex* FindSimpleIP(AliVEvent* event,
944  Double_t maxDispersion=0.04,
945  Double_t maxZError=0.25);
955  const AliVVertex* FindRealIP(AliVEvent* event,
956  Double_t maxDispersion=0.04,
957  Double_t maxZError=0.25);
967  const AliVVertex* CheckIP(const AliVVertex* ip,
968  Double_t maxDispersion=0.04,
969  Double_t maxZError=0.25);
970 
980  const AliVVertex* FindIP(AliVEvent* event,
981  Double_t maxDispersion=0.04,
982  Double_t maxZError=0.25);
993  Double_t FindMultCentrality(AliVEvent* event, Int_t& nTracklets);
1003  Double_t FindCompatCentrality(AliVEvent* event);
1015  Double_t FindCentrality(AliVEvent* event, Int_t& nTracklets);
1021  Bool_t FindTrigger();
1022  /* @} */
1023 
1024  // -----------------------------------------------------------------
1038  virtual Double_t LookupWeight(AliAODTracklet* tracklet,
1039  Double_t cent,
1040  Double_t ipz);
1050  virtual UShort_t CheckTracklet(AliAODTracklet* tracklet) const;
1051  /* @} */
1052 
1053  // -----------------------------------------------------------------
1065  void ProcessEvent(Double_t cent,const AliVVertex* ip,TClonesArray* tracklets);
1066  /* @} */
1067 
1068  // -----------------------------------------------------------------
1080  virtual Bool_t MasterFinalize(Container* results);
1081  /* @} */
1082 
1083  // -----------------------------------------------------------------
1085  enum {
1086  kAll=1, // Count all events
1087  kEvent, // Have event
1088  kTracklets, // Have tracklets
1089  kTrigger, // Have trigger
1090  kIP, // Have IP
1091  kCentrality, // Have centrality
1092  kCompleted // Have completed
1093  };
1094 
1095  // -----------------------------------------------------------------
1102 
1104 
1106 
1108 
1109  TProfile* fNBareVsGood;
1111  TProfile* fNBareVsFake;
1113  TProfile* fNTrackletVsGood;
1115  TProfile* fNTrackletVsFake;
1119  TProfile* fNGeneratedVsFake;
1121  TProfile* fCentTracklets;
1122 
1123  TProfile* fCentEst;
1124 
1152 
1154 };
1155 //====================================================================
1162 {
1163 public:
1165  typedef TList Container;
1171  {}
1175  AliTrackletAODMCdNdeta(const char* name)
1176  : AliTrackletAODdNdeta(name)
1177  {}
1185  {}
1190  //__________________________________________________________________
1207  {
1213  fCombinatorics(0),
1214  fDistinct(0),
1215  fPrimaries(0),
1216  fSecondaries(0),
1217  fGenerated(0)
1218  {
1219  }
1226  CentBin(Double_t c1, Double_t c2);
1232  CentBin(const CentBin& o)
1234  fCombinatorics(0),
1235  fPrimaries(0),
1236  fSecondaries(0),
1237  fGenerated(0)
1238  {}
1242  virtual ~CentBin() {}
1248  CentBin& operator=(const CentBin&) { return *this; }
1249  protected:
1255 
1256  ClassDef(CentBin,1);
1257  };
1258 
1259 protected:
1260  // -----------------------------------------------------------------
1274  {
1275  return new CentBin(c1, c2);
1276  }
1277  /* @} */
1278 
1279  // -----------------------------------------------------------------
1284  virtual const char* GetBranchName() const { return "AliAODMCTracklets"; }
1285  /* @} */
1286 
1288 };
1289 //====================================================================
1296 {
1297 public:
1299  typedef TList Container;
1305  fWeights(0),
1306  fEtaWeight(0),
1307  fWeightCorr(0)
1308  {}
1313  : AliTrackletAODMCdNdeta(name),
1314  fWeights(0),
1315  fEtaWeight(0),
1316  fWeightCorr(0)
1317  {}
1325  fWeights(0),
1326  fEtaWeight(0),
1327  fWeightCorr(0)
1328  {}
1347  void Print(Option_t* option="") const;
1348 
1349  // -----------------------------------------------------------------
1365  void SetWeightCalc(UChar_t mode=0) {
1366  Info("SetWeightCalc", "Use=%d", mode);
1367  if (fWeights) fWeights->SetCalc(mode);}
1373  void SetWeightMask(UChar_t mask=0xFF) {
1374  Info("SetWeightMask", "mask=0x%x", mask);
1375  if (fWeights) fWeights->SetMask(mask); }
1381  void SetWeightVeto(UChar_t veto=0x0) {
1382  Info("SetWeightVeto", "veto=0x%x", veto);
1383  if (fWeights) fWeights->SetVeto(veto); }
1392  /* @} */
1393 protected:
1394  // -----------------------------------------------------------------
1404  Bool_t WorkerInit();
1414  virtual Double_t LookupWeight(AliAODTracklet* tracklet,
1415  Double_t cent,
1416  Double_t ipz);
1417  /* @} */
1418  // -----------------------------------------------------------------
1430  Bool_t MasterFinalize(Container* results);
1431  /* @} */
1432  // -----------------------------------------------------------------
1434  TProfile2D* fEtaWeight;
1437 };
1438 
1439 //____________________________________________________________________
1441  : AliAnalysisTaskSE(),
1442  fContainer(0),
1443  fCentBins(0),
1444  fIPz(0),
1445  fCent(0),
1446  fStatus(0),
1447  fEtaPhi(0),
1448  fNBareVsGood(0),
1449  fNBareVsFake(0),
1450  fNTrackletVsGood(0),
1451  fNTrackletVsFake(0),
1452  fNGeneratedVsGood(0),
1453  fNGeneratedVsFake(0),
1454  fCentTracklets(0),
1455  fCentEst(0),
1456  fCentMethod(""),
1457  fCentIdx(-1),
1458  fTrkIdx(-1),
1459  fCentAxis(1,0,0),
1460  fIPzAxis(1,0,0),
1461  fEtaAxis(1,0,0),
1462  fPhiAxis(1,0,0),
1463  fMaxDelta(0),
1464  fTailDelta(0),
1465  fTailMax(-1),
1466  fDPhiShift(0),
1467  fShiftedDPhiCut(0),
1468  fDeltaCut(0),
1469  fAbsMinCent(-1)
1470 {}
1471 //____________________________________________________________________
1473  : AliAnalysisTaskSE(name),
1474  fContainer(0),
1475  fCentBins(0),
1476  fIPz(0),
1477  fCent(0),
1478  fStatus(0),
1479  fEtaPhi(0),
1480  fNBareVsGood(0),
1481  fNBareVsFake(0),
1482  fNTrackletVsGood(0),
1483  fNTrackletVsFake(0),
1484  fNGeneratedVsGood(0),
1485  fNGeneratedVsFake(0),
1486  fCentTracklets(0),
1487  fCentEst(0),
1488  fCentMethod("V0M"),
1489  fCentIdx(-1),
1490  fTrkIdx(-1),
1491  fCentAxis(10,0,100),
1492  fIPzAxis(30,-15,+15),
1493  fEtaAxis(16,-2,+2),
1494  fPhiAxis(100,0,TMath::TwoPi()),
1495  fMaxDelta(25),
1496  fTailDelta(5),
1497  fTailMax(-1),
1498  fDPhiShift(0.0045),
1499  fShiftedDPhiCut(-1),
1500  fDeltaCut(1.5),
1501  fAbsMinCent(-1)
1502 {
1503  FixAxis(fCentAxis, "Centrality [%]");
1504  FixAxis(fIPzAxis, "IP_{#it{z}} [cm]");
1505  FixAxis(fEtaAxis, "#eta");
1506  FixAxis(fPhiAxis, "#phi");
1507 
1508  DefineOutput(1, Container::Class());
1509  DefineOutput(2, Container::Class());
1510 }
1511 //____________________________________________________________________
1513  : AliAnalysisTaskSE(o),
1514  fContainer(0),
1515  fCentBins(0),
1516  fIPz(0),
1517  fCent(0),
1518  fStatus(0),
1519  fEtaPhi(0),
1520  fNBareVsGood(0),
1521  fNBareVsFake(0),
1522  fNTrackletVsGood(0),
1523  fNTrackletVsFake(0),
1524  fNGeneratedVsGood(0),
1525  fNGeneratedVsFake(0),
1526  fCentTracklets(0),
1527  fCentEst(0),
1528  fCentMethod(o.fCentMethod),
1529  fCentIdx(o.fCentIdx),
1530  fTrkIdx(o.fTrkIdx),
1531  fCentAxis(o.fCentAxis),
1532  fIPzAxis(o.fIPzAxis),
1533  fEtaAxis(o.fEtaAxis),
1534  fPhiAxis(o.fPhiAxis),
1535  fMaxDelta(o.fMaxDelta),
1536  fTailDelta(o.fTailDelta),
1537  fTailMax(o.fTailMax),
1538  fDPhiShift(o.fDPhiShift),
1539  fShiftedDPhiCut(o.fShiftedDPhiCut),
1540  fDeltaCut(o.fDeltaCut),
1541  fAbsMinCent(o.fAbsMinCent)
1542 {}
1543 //____________________________________________________________________
1546 {
1547  if (&o == this) return *this;
1548  if (fContainer) {
1549  delete fContainer;
1550  fContainer = 0;
1551  }
1552  if (fCentBins) {
1553  delete fCentBins;
1554  fCentBins = 0;
1555  }
1556  fIPz = 0;
1557  fCent = 0;
1559  fCentIdx = o.fCentIdx;
1560  fTrkIdx = o.fTrkIdx;
1561  fCentAxis = o.fCentAxis;
1562  fIPzAxis = o.fIPzAxis;
1563  fEtaAxis = o.fEtaAxis;
1564  fPhiAxis = o.fPhiAxis;
1565  fMaxDelta = o.fMaxDelta;
1566  fTailDelta = o.fTailDelta;
1567  fTailMax = o.fTailMax;
1568  fDPhiShift = o.fDPhiShift;
1570  fDeltaCut = o.fDeltaCut;
1572  return *this;
1573 }
1574 //____________________________________________________________________
1577  o)
1578 {
1579  if (&o == this) return *this;
1581  return *this;
1582 }
1583 //____________________________________________________________________
1585  const char* resFile)
1586 {
1587  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1588  if (!mgr) {
1589  AliError("No analysis manager to connect to.");
1590  return false;
1591  }
1592 
1593  // Add to the manager
1594  mgr->AddTask(this);
1595 
1596  // Create and connect output containers
1597  TString sumOut;
1598  TString resOut;
1599  if (sumFile && sumFile[0] != '\0') sumOut = sumFile;
1600  if (resFile && resFile[0] != '\0') resOut = resFile;
1601  else if (sumFile && sumFile[0] != '\0') resOut = sumFile;
1602 
1603  // If the string is null or 'default' connect to standard output file
1604  if (sumOut.IsNull() || sumOut.EqualTo("default", TString::kIgnoreCase))
1605  sumOut = AliAnalysisManager::GetCommonFileName();
1606  // If the string is null or 'default' connect to standard output file
1607  if (resOut.IsNull() || resOut.EqualTo("default", TString::kIgnoreCase))
1608  resOut = AliAnalysisManager::GetCommonFileName();
1609 
1610  // Always connect input
1611  mgr->ConnectInput(this, 0, mgr->GetCommonInputContainer());
1612 
1613  // Connect sum list unless the output 'none' is specified
1614  if (!sumOut.EqualTo("none", TString::kIgnoreCase)) {
1615  TString sumName(Form("%sSums", GetName()));
1616  AliAnalysisDataContainer* sumCon =
1617  mgr->CreateContainer(sumName, TList::Class(),
1618  AliAnalysisManager::kOutputContainer, sumOut);
1619  mgr->ConnectOutput(this, 1, sumCon);
1620  }
1621 
1622  // Connect the result list unless the output 'none' is specified
1623  if (!resOut.EqualTo("none", TString::kIgnoreCase)) {
1624  TString resName(Form("%sResults", GetName()));
1625  AliAnalysisDataContainer* resCon =
1626  mgr->CreateContainer(resName, TList::Class(),
1627  AliAnalysisManager::kParamContainer, resOut);
1628  mgr->ConnectOutput(this, 2, resCon);
1629  }
1630  return true;
1631 }
1632 //____________________________________________________________________
1634 {
1635  Double_t shiftedDPhiCut = fShiftedDPhiCut;
1636  if (shiftedDPhiCut < 0) shiftedDPhiCut = TMath::Sqrt(fDeltaCut)*0.06;
1637  Double_t tailMax = fTailMax;
1638  if (tailMax < 0) tailMax = fMaxDelta;
1639 
1640  Printf("%s: %s", ClassName(), GetName());
1641  Printf(" %22s: 0x%08x", "Off-line trigger mask", fOfflineTriggerMask);
1642  Printf(" %22s: %f", "Delta phi shift", fDPhiShift);
1643  Printf(" %22s: %f", "Shifted Delta phi cut", shiftedDPhiCut);
1644  Printf(" %22s: %f", "Delta cut", fDeltaCut);
1645  Printf(" %22s: %f", "max Delta", fMaxDelta);
1646  Printf(" %22s: %f", "tail Delta", fTailDelta);
1647  Printf(" %22s: %f", "tail maximum", tailMax);
1648  Printf(" %22s: %f%%", "Absolute least c", fAbsMinCent);
1651  PrintAxis(fIPzAxis,1,"IPz");
1652  PrintAxis(fCentAxis,0);
1653 
1654  if (!fCentBins) return;
1655 
1656  Printf("--- Centrality bins");
1657  TIter next(fCentBins);
1658  CentBin* bin = 0;
1659  while ((bin = static_cast<CentBin*>(next()))) {
1660  bin->Print(option);
1661  }
1662 }
1663 //____________________________________________________________________
1665 {
1667  if (!fWeights) return;
1668  fWeights->Print(option);
1669 }
1670 
1671 //____________________________________________________________________
1673 {
1674  Printf(" Centrality bin: %s", fName.Data());
1675  Printf(" Low cut: %5.1f", fLow);
1676  Printf(" High cut: %5.1f", fHigh);
1677 
1678  if (!fSubs) return;
1679  Printf(" --- Histogram sets");
1680  TIter next(fSubs);
1681  Histos* h = 0;
1682  while ((h = static_cast<Histos*>(next()))) {
1683  h->Print(option);
1684  }
1685 }
1686 
1687 //____________________________________________________________________
1689 {
1690  Sub::SetDebug(lvl);
1691  TIter next(fSubs);
1692  Histos* h = 0;
1693  while ((h = static_cast<Histos*>(next()))) {
1694  h->SetDebug(lvl);
1695  }
1696 }
1697 namespace {
1698  void Bits2String(UChar_t m, char out[7])
1699  {
1700  if (m & AliAODTracklet::kInjection) out[0] = 'I'; else out[0] = '-';
1701  if (m & AliAODTracklet::kCombinatorics) out[1] = 'C'; else out[1] = '-';
1702  if (m & AliAODTracklet::kSecondary) out[2] = 'S'; else out[2] = '-';
1703  if (m & AliAODTracklet::kDistinct) out[3] = 'D'; else out[3] = '-';
1704  if (m & AliAODTracklet::kSimulated) out[4] = 'X'; else out[4] = '-';
1705  if (m & AliAODTracklet::kGenerated) out[5] = 'G'; else out[5] = '-';
1706  out[6] = '\0';
1707  }
1708 }
1709 
1710 //____________________________________________________________________
1712 {
1713  char cMask[7]; Bits2String(fMask, cMask);
1714  char cVeto[7]; Bits2String(fVeto, cVeto);
1715  Printf(" Histograms: %s", fName.Data());
1716  Printf(" Mask: 0x%02x (%s)", fMask, cMask);
1717  Printf(" Veto: 0x%02x (%s)", fVeto, cVeto);
1718  Printf(" Delta: %s", fEtaDeltaIPz ? "yes" : "no");
1719  Printf(" Delta per PDG: %s", fEtaDeltaPdg ? "yes" : "no");
1720 }
1721 
1722 
1723 //____________________________________________________________________
1724 void
1726 {
1727  if (!WorkerInit()) {
1728  AliWarning("Failed to initialize on worker");
1729  return;
1730  }
1731  PostData(1,fContainer);
1732 }
1733 
1734 //____________________________________________________________________
1736 {
1737  if (DebugLevel() > 1) Printf("Initialising on worker");
1738  if (fShiftedDPhiCut < 0) fShiftedDPhiCut = TMath::Sqrt(fDeltaCut)*0.06;
1739  if (fTailMax < 0) fTailMax = fMaxDelta;
1740 
1741  fContainer = new Container;
1742  fContainer->SetName(Form("%sSums", GetName()));
1743  fContainer->SetOwner();
1744 
1745  Bool_t save = TH1::GetDefaultSumw2();
1746  TH1::SetDefaultSumw2();
1747  fIPz = Make1D(fContainer, "ipz", "", kMagenta+2, 20, fIPzAxis);
1748  fCent = Make1D(fContainer, "cent", "", kMagenta+2, 20, fCentAxis);
1749  fEtaPhi = Make2D(fContainer, "etaPhi","",kMagenta+2, 20, fEtaAxis,fPhiAxis);
1750 
1751  TAxis trackletAxis(1000, 0, 10000);
1752  FixAxis(trackletAxis, "#it{N}_{tracklets}");
1753  fNBareVsGood = Make1P(fContainer, "nBareVsGood", "Good",
1754  kGreen+2, 20, trackletAxis);
1755  fNBareVsFake = Make1P(fContainer, "nBareVsFake", "Fake",
1756  kRed+2, 20, trackletAxis);
1757  fNBareVsGood->SetYTitle("#LT#it{N}_{tracklet}#GT");
1758  fNBareVsFake->SetYTitle("#LT#it{N}_{tracklet}#GT");
1759  fNBareVsGood->SetStats(0);
1760  fNBareVsFake->SetStats(0);
1761 
1762  fNTrackletVsGood = Make1P(fContainer, "nTrackletVsGood", "Good",
1763  kGreen+2, 20, trackletAxis);
1764  fNTrackletVsFake = Make1P(fContainer, "nTrackletVsFake", "Fake",
1765  kRed+2, 20, trackletAxis);
1766  fNTrackletVsGood->SetYTitle("#LT#it{N}_{tracklet}#GT");
1767  fNTrackletVsFake->SetYTitle("#LT#it{N}_{tracklet}#GT");
1768  fNTrackletVsGood->SetStats(0);
1769  fNTrackletVsFake->SetStats(0);
1770 
1771  TAxis generatedAxis(1000, 0, 15000);
1772  FixAxis(generatedAxis, "#it{N}_{generated,|#eta|<2}");
1773  fNGeneratedVsGood = Make1P(fContainer, "nGeneratedVsGood", "Good",
1774  kGreen+2, 24, generatedAxis);
1775  fNGeneratedVsFake = Make1P(fContainer, "nGeneratedVsFake", "Fake",
1776  kRed+2, 24, generatedAxis);
1777  fNGeneratedVsGood->SetYTitle("#LT#it{N}_{tracklet}#GT");
1778  fNGeneratedVsFake->SetYTitle("#LT#it{N}_{tracklet}#GT");
1779  fNGeneratedVsGood->SetStats(0);
1780  fNGeneratedVsFake->SetStats(0);
1781 
1782  fCentTracklets = new TProfile("centTracklets",
1783  "Mean number of tracklets per centrality",
1784  1050, 0, 105);
1785  fCentTracklets->SetXTitle("Centrality [%]");
1786  fCentTracklets->SetYTitle("#LTtracklets#GT");
1787  fCentTracklets->SetDirectory(0);
1788  fCentTracklets->SetMarkerStyle(20);
1789  fCentTracklets->SetMarkerColor(kMagenta+2);
1790  fCentTracklets->SetLineColor(kMagenta+2);
1791  fContainer->Add(fCentTracklets);
1792  fCentEst = Make1P(fContainer,"centEstimator","",kMagenta+2, 20, fCentAxis);
1793  fCentEst->SetYTitle(Form("#LT%s%GT",fCentMethod.Data()));
1794  TH1::SetDefaultSumw2(save);
1795 
1796  fStatus = new TH1F("status", "Status of task",
1797  kCompleted, .5, kCompleted+.5);
1798  fStatus->SetMarkerSize(2);
1799  fStatus->SetMarkerColor(kMagenta+2);
1800  fStatus->SetLineColor(kMagenta+2);
1801  fStatus->SetFillColor(kMagenta+2);
1802  fStatus->SetFillStyle(1001);
1803  fStatus->SetBarOffset(0.1);
1804  fStatus->SetBarWidth(0.4);
1805  fStatus->SetDirectory(0);
1806  fStatus->SetStats(0);
1807  fStatus->SetXTitle("Event have");
1808  fStatus->SetYTitle("# Events");
1809  fStatus->GetXaxis()->SetBinLabel(kAll, "Been seen");
1810  fStatus->GetXaxis()->SetBinLabel(kEvent, "Event data");
1811  fStatus->GetXaxis()->SetBinLabel(kTracklets, "Tracklets");
1812  fStatus->GetXaxis()->SetBinLabel(kTrigger, "Trigger");
1813  fStatus->GetXaxis()->SetBinLabel(kIP, "IP");
1814  fStatus->GetXaxis()->SetBinLabel(kCentrality, "Centrality");
1815  fStatus->GetXaxis()->SetBinLabel(kCompleted, "Completed");
1816  fContainer->Add(fStatus);
1817 
1818  typedef TParameter<double> DP;
1819  typedef TParameter<bool> BP;
1820  typedef TParameter<int> IP;
1821  Container* params = new Container;
1822  params->SetName("parameters");
1823  params->SetOwner();
1824  fContainer->Add(params);
1825  params->Add(new DP("DPhiShift", fDPhiShift, 'f'));
1826  params->Add(new DP("ShiftedDPhiCut", fShiftedDPhiCut, 'f'));
1827  params->Add(new DP("DeltaCut", fDeltaCut, 'f'));
1828  params->Add(new DP("MaxDelta", fMaxDelta, 'f'));
1829  params->Add(new DP("TailDelta", fTailDelta, 'f'));
1830  params->Add(new DP("TailMax", fTailMax, 'f'));
1831  params->Add(new DP("AbsMinCent", fAbsMinCent, 'f'));
1832  // Create our centrality bins
1833  if (!InitCentBins(0)) {
1834  AliWarning("Failed to initialize centrality bins");
1835  return false;
1836  }
1837 
1838  // Print information to log
1839  Print();
1840  return true;
1841 }
1842 //____________________________________________________________________
1844 {
1845  if (!AliTrackletAODMCdNdeta::WorkerInit()) return false;
1846  fEtaWeight = 0;
1847  if (!fWeights) {
1848  AliFatal("No weights set!");
1849  return false;
1850  }
1851 
1852  TAxis wAxis(100,0,10);
1853  FixAxis(wAxis,"#it{w}_{i}");
1854  fEtaWeight = Make2P(fContainer, "etaWeight", "#LTw#GT", kYellow+2, 24,
1855  fEtaAxis, fCentAxis);
1856  fWeightCorr = Make2D(fContainer, "weightCorr", "w_{1} vs w_{2}",
1857  kCyan+2, 24, wAxis, wAxis);
1858  fWeights->Store(fContainer);
1859  return true;
1860 }
1861 //____________________________________________________________________
1864  Float_t c1,
1865  Float_t c2,
1866  Container* existing,
1867  const TAxis& deltaAxis)
1868 {
1869  CentBin* ret = MakeCentBin(c1, c2);
1870  Bool_t ok = true;
1871  if (!existing)
1872  ok = ret->WorkerInit(fContainer,fEtaAxis,fIPzAxis,deltaAxis);
1873  else
1874  ok = ret->FinalizeInit(existing);
1875  if (!ok) {
1876  AliWarningF("Failed to initialize bin %s", ret->GetName());
1877  delete ret;
1878  return 0;
1879  }
1880  ret->SetDebug(DebugLevel());
1881  fCentBins->AddAt(ret, bin);
1882  return ret;
1883 
1884 }
1885 //____________________________________________________________________
1887 {
1888  if (DebugLevel() > 1)
1889  Printf("Initialising on centrality bins on %s",
1890  existing ? "master" : "worker");
1891  if (fCentBins) return true;
1892  fCentBins = new Container;
1893  fCentBins->SetName("centralityBins");
1894  fCentBins->SetOwner();
1895 
1896  TAxis deltaAxis (Int_t(5*fMaxDelta),0,fMaxDelta);
1897  FixAxis(deltaAxis,
1898  "#Delta=[(#Delta#phi-#delta#phi)/#sigma_{#phi}]^{2}+"
1899  "[#Delta#thetasin^{-2}(#theta)/#sigma_{#theta}]^{2}");
1900 
1901  // Add min-bias bin
1902  if (!InitCentBin(0,0,0,existing,deltaAxis)) return false;
1903 
1904  // Add other bins
1905  Int_t nCentBins = fCentAxis.GetNbins();
1906  for (Int_t i = 1; i <= nCentBins; i++) {
1907  Float_t c1 = fCentAxis.GetBinLowEdge(i);
1908  Float_t c2 = fCentAxis.GetBinUpEdge(i);
1909  if (!InitCentBin(i, c1, c2, existing, deltaAxis)) return false;
1910  }
1911  if (!InitCentBin(nCentBins+1, 0, 100, existing, deltaAxis)) return false;
1912  return true;
1913 }
1914 
1915 //____________________________________________________________________
1917  : Sub(""),
1918  fSubs(0),
1919  fLow(c1),
1920  fHigh(c2),
1921  fIPz(0),
1922  fCent(0),
1923  fCentIPz(0),
1924  fMeasured(0),
1925  fInjection(0)
1926 {
1927  if (c1 >= c2)
1928  fName = "all";
1929  else
1931  fMeasured = MakeHistos("measured", kRed+2, 20,
1932  kMeasuredMask, // No requirements, just veto
1933  kMeasuredVeto);
1934  fInjection = MakeHistos("injected", kOrange+2, 21,
1935  kInjectedMask,
1936  kInjectedVeto);
1937  fSubs = new Container;
1938  fSubs->SetOwner(true);
1939  fSubs->Add(fMeasured);
1940  fSubs->Add(fInjection);
1941 }
1942 
1943 //____________________________________________________________________
1945  : AliTrackletAODdNdeta::CentBin(c1, c2),
1946  fCombinatorics(0),
1947  fDistinct(0),
1948  fPrimaries(0),
1949  fSecondaries(0),
1950  fGenerated(0)
1951 {
1952  // Combinatorics is everything that has the combinatorics bit
1953  // Primaries all with simulated bit, but not secondary or
1954  // combinatorics bit.
1955  // Secondaries are all those with the secondary bit set
1956  fCombinatorics = MakeHistos("combinatorics", kMagenta+2, 30,
1959  fDistinct = MakeHistos("distinct", kCyan+2, 27,
1960  kDistinctMask,
1961  kDistinctVeto);
1962  fPrimaries = MakeHistos("primaries", kGreen+2, 26,
1963  kPrimaryMask,
1964  kPrimaryVeto);
1965  fSecondaries = MakeHistos("secondaries", kBlue+2, 32,
1967  kSecondaryVeto);
1968  fGenerated = MakeHistos("generated", kGray+1, 28,
1970  kGeneratedVeto);
1971  fMeasured->fStyle = 24;
1972  fInjection->fStyle = 25;
1973  fSubs->Add(fCombinatorics);
1974  fSubs->Add(fDistinct);
1975  fSubs->Add(fPrimaries);
1976  fSubs->Add(fSecondaries);
1977  fSubs->AddAfter(fMeasured, fGenerated);
1978 }
1979 
1980 //____________________________________________________________________
1983  Color_t color,
1984  Style_t style,
1985  UShort_t mask,
1986  UShort_t veto)
1987 {
1988  return new Histos(name, color, style, mask, veto);
1989 }
1990 
1991 //____________________________________________________________________
1993  const TAxis& etaAxis,
1994  const TAxis& ipzAxis,
1995  const TAxis& deltaAxis)
1996 {
1997  if (fDebug > 0)
1998  Printf("Initializing centrality bin %s", fName.Data());
1999  if (!Sub::WorkerInit(parent, etaAxis, ipzAxis, deltaAxis)) return false;
2000 
2001  TAxis centAxis(20, fLow, fHigh);
2002  FixAxis(centAxis, "Centrality [%]");
2003  fCent = Make1D(fContainer,"cent","Centrality [%]",
2004  kMagenta+2,20,centAxis);
2005  fIPz = Make1D(fContainer,"ipz","IP_{#it{z}} [cm]",kRed+2,20,ipzAxis);
2006  fCentIPz = Make1P(fContainer,"centIpz","#LTc#GT vs IP_{#it{z}}",kPink+2,
2007  20,ipzAxis);
2008  fCentIPz->SetYTitle("#LTc#GT");
2009  TIter next(fSubs);
2010  Histos* h = 0;
2011  while ((h = static_cast<Histos*>(next()))) {
2012  if (!h ->WorkerInit(fContainer, etaAxis,ipzAxis,deltaAxis))
2013  return false;
2014  }
2015 
2016  return true;
2017 }
2018 //____________________________________________________________________
2020  const TAxis& etaAxis,
2021  const TAxis& ipzAxis,
2022  const TAxis& deltaAxis)
2023 {
2024  if (!Sub::WorkerInit(parent, etaAxis, ipzAxis, deltaAxis)) return false;
2025  TString shrt(Form("%c", GetName()[0]));
2026  shrt.ToUpper();
2027  if (GetC(parent, "generated", false) != 0) shrt.Append("'");
2028 
2029  // Do not make eta vs IPz for secondaries and primaries
2030  if (!IsPrimary() && !IsSecondary())
2031  fEtaIPz = Make2D(fContainer, "etaIPz", shrt.Data(),
2032  kRed+2, 20, etaAxis, ipzAxis);
2033  // Always make eta vs Delta distribution, except for MC truth
2034  if (!IsGenerated())
2035  fEtaDeltaIPz = Make3D(fContainer, "etaDeltaIPz",
2036  Form("#Delta_{%s}",shrt.Data()),
2037  kBlue+2, 21, etaAxis, deltaAxis, ipzAxis);
2038  if (!IsGenerated() &&
2039  // !IsPrimary() &&
2040  // !IsSecondary() &&
2041  !IsInjected() &&
2042  fStyle != 20) // Last condition to not make this for real data
2043  fEtaPdgIPz = Make3D(fContainer, "etaPdgIPz",
2044  "Parent particle type",
2045  kGreen+2, 22, etaAxis, PdgAxis(), ipzAxis);
2046 
2047  if (IsGenerated()) {
2048  fEtaPdg = Make2D(fContainer, "etaPdg",
2049  "Primary particle type",
2050  kYellow+2, 30, etaAxis, PdgAxis());
2051  TAxis ptAxis(100, 0, 5);
2052  ptAxis.SetTitle("#it{p}_{T}");
2053  FixAxis(ptAxis);
2054  fEtaPt = Make2D(fContainer, "etaPt",
2055  "Primary transverse momentum",
2056  kCyan+2, 30, etaAxis, ptAxis);
2057  }
2058  if (IsPrimary() || IsSecondary() || IsCombinatoric()) {
2059  fEtaDeltaPdg = Make3D(fContainer, "etaDeltaPdg",
2060  "#Delta by primary particle type",
2061  kMagenta, 22, etaAxis, deltaAxis, PdgAxis());
2062  }
2063 
2064  SetAttr(fColor, fStyle);
2065  return true;
2066 }
2067 //____________________________________________________________________
2069 {
2070  if (fEtaIPz) {
2071  fEtaIPz->SetMarkerStyle(m);
2072  fEtaIPz->SetMarkerColor(c);
2073  fEtaIPz->SetLineColor(c);
2074  fEtaIPz->SetFillColor(c);
2075  }
2076  if (fEtaDeltaIPz) {
2077  fEtaDeltaIPz->SetMarkerStyle(m);
2078  fEtaDeltaIPz->SetMarkerColor(c);
2079  fEtaDeltaIPz->SetLineColor(c);
2080  fEtaDeltaIPz->SetFillColor(c);
2081  }
2082  if (fEtaDeltaPdg) {
2083  fEtaDeltaPdg->SetMarkerStyle(m);
2084  fEtaDeltaPdg->SetMarkerColor(c);
2085  fEtaDeltaPdg->SetLineColor(c);
2086  fEtaDeltaPdg->SetFillColor(c);
2087  }
2088  if (fEtaPdgIPz) {
2089  fEtaPdgIPz->SetMarkerStyle(m);
2090  fEtaPdgIPz->SetMarkerColor(c);
2091  fEtaPdgIPz->SetLineColor(c);
2092  fEtaPdgIPz->SetFillColor(c);
2093  }
2094  if (fEtaPdg) {
2095  fEtaPdg->SetMarkerStyle(m);
2096  fEtaPdg->SetMarkerColor(c);
2097  fEtaPdg->SetLineColor(c);
2098  fEtaPdg->SetFillColor(c);
2099  }
2100  if (fEtaPt) {
2101  fEtaPt->SetMarkerStyle(m);
2102  fEtaPt->SetMarkerColor(c);
2103  fEtaPt->SetLineColor(c);
2104  fEtaPt->SetFillColor(c);
2105  }
2106 }
2107 
2108 //____________________________________________________________________
2109 void
2111 {
2112  if (DebugLevel() > 0) Printf("In user exec");
2113  Double_t cent = -1;
2114  const AliVVertex* ip = 0;
2115  TClonesArray* tracklets = 0;
2116  if (!CheckEvent(cent, ip, tracklets)) {
2117  AliWarningF("Event didn't pass %f, %p, %p", cent, ip, tracklets);
2118  Printf("Argh, check data failed %f, %p, %p", cent, ip, tracklets);
2119  return;
2120  }
2121  if (DebugLevel() > 0) Printf("Got centrality=%f ipZ=%f %d tracklets",
2122  cent, ip->GetZ(), tracklets->GetEntriesFast());
2123  ProcessEvent(cent, ip, tracklets);
2124 
2125  PostData(1,fContainer);
2126  fStatus->Fill(kCompleted);
2127 }
2128 
2129 //____________________________________________________________________
2131  const AliVVertex*& ip,
2132  TClonesArray*& tracklets)
2133 {
2134  // Count all events
2135  fStatus->Fill(kAll);
2136 
2137  // Check for event
2138  AliVEvent* event = InputEvent();
2139  if (!event) {
2140  AliWarning("No event");
2141  return false;
2142  }
2143  fStatus->Fill(kEvent);
2144 
2145  // Check if we have the tracklets
2146  tracklets = FindTracklets(event);
2147  if (!tracklets) return false;
2148  fStatus->Fill(kTracklets);
2149 
2150  // Check if event was triggered
2151  Bool_t trg = FindTrigger();
2152  if (!trg) return false;
2153  fStatus->Fill(kTrigger);
2154 
2155  // Check the interaction point
2156  ip = FindIP(event);
2157  if (!ip) return false;
2158  fStatus->Fill(kIP);
2159 
2160  // Check the centrality
2161  Int_t nTracklets = 0;
2162  cent = FindCentrality(event, nTracklets);
2163  // Do not fail on missing centrality
2164  if (cent >= 0) fStatus->Fill(kCentrality);
2165 
2166  fIPz->Fill(ip->GetZ());
2167  fCent->Fill(cent);
2168  fCentTracklets->Fill(cent, nTracklets);
2169 
2170  return true;
2171 }
2172 
2173 //____________________________________________________________________
2174 TClonesArray* AliTrackletAODdNdeta::FindTracklets(AliVEvent* event)
2175 {
2176  // Check the multiplicity
2177  TObject* obj = event->FindListObject(GetBranchName());
2178  if (!obj) {
2179  AliWarningF("Couldn't get object %s", GetBranchName());
2180  // event->GetList()->Print();
2181  return 0;
2182  }
2183  if (!obj->IsA()->InheritsFrom(TClonesArray::Class())) {
2184  AliWarningF("Object %s is not a TClonesArray but a %s",
2185  obj->GetName(), obj->ClassName());
2186  return 0;
2187  }
2188  return static_cast<TClonesArray*>(obj);
2189 }
2190 
2191 //____________________________________________________________________
2192 const AliVVertex* AliTrackletAODdNdeta::CheckIP(const AliVVertex* ip,
2193  Double_t maxDispersion,
2194  Double_t maxZError)
2195 {
2196  if (!ip) return 0;
2197  if (ip->GetNContributors() <= 0) {
2198  AliWarning("Not enough contributors for IP");
2199  return 0;
2200  }
2201  // If this is from the Z vertexer, do some checks
2202  if (ip->IsFromVertexerZ()) {
2203  // Get covariance matrix
2204  Double_t covar[6];
2205  ip->GetCovarianceMatrix(covar);
2206  Double_t sigmaZ = TMath::Sqrt(covar[5]);
2207  if (sigmaZ >= maxZError) {
2208  AliWarningF("IPz resolution = %f >= %f", sigmaZ, maxZError);
2209  return 0;
2210  }
2211 
2212  // If this IP doesn not derive from AliVertex, don't check dispersion.
2213  if (ip->IsA()->InheritsFrom(AliVertex::Class())) {
2214  const AliVertex* ipv = static_cast<const AliVertex*>(ip);
2215  // Dispersion is the parameter used by the vertexer for finding the IP.
2216  if (ipv->GetDispersion() >= maxDispersion) {
2217  AliWarningF("IP dispersion = %f >= %f",
2218  ipv->GetDispersion(), maxDispersion);
2219  return 0;
2220  }
2221  }
2222  }
2223 
2224  // If we get here, we either have a full 3D vertex or track
2225  // vertex, and we should check if it is in range
2226  if (ip->GetZ() < fIPzAxis.GetXmin() || ip->GetZ() > fIPzAxis.GetXmax()) {
2227  AliWarningF("IPz = %fcm out of range [%f,%f]cm",
2228  ip->GetZ(), fIPzAxis.GetXmin(), fIPzAxis.GetXmax());
2229  return 0;
2230  }
2231  return ip;
2232 }
2233 
2234 //____________________________________________________________________
2235 const AliVVertex* AliTrackletAODdNdeta::FindRealIP(AliVEvent* event,
2236  Double_t maxDispersion,
2237  Double_t maxZError)
2238 {
2239  const AliVVertex* ip = event->GetPrimaryVertex();
2240  if (!ip) {
2241  if (DebugLevel() > 1) AliWarning("No real IP for this event found!");
2242  return 0;
2243  }
2244  return CheckIP(ip, maxDispersion, maxZError);
2245 }
2246 
2247 //____________________________________________________________________
2248 const AliVVertex* AliTrackletAODdNdeta::FindSimpleIP(AliVEvent* event,
2249  Double_t maxDispersion,
2250  Double_t maxZError)
2251 {
2252  static AliVertex* ip;
2253  AliAODSimpleHeader* head =
2254  static_cast<AliAODSimpleHeader*>(event
2255  ->FindListObject("AliAODSimpleHeader"));
2256  if (!head) {
2257  if (DebugLevel() > 1) AliWarning("No simple header");
2258  return 0;
2259  }
2260  if (!ip)
2261  ip = new AliVertex;
2262  ip->SetXv(head->fRecIP.X());
2263  ip->SetYv(head->fRecIP.Y());
2264  ip->SetZv(head->fRecIP.Z());
2265  ip->SetNContributors(10);
2266  ip->SetTitle("");
2267  return CheckIP(ip, maxDispersion, maxZError);
2268 }
2269 
2270 //____________________________________________________________________
2271 const AliVVertex* AliTrackletAODdNdeta::FindIP(AliVEvent* event,
2272  Double_t maxDispersion,
2273  Double_t maxZError)
2274 {
2275  const AliVVertex* ip = FindRealIP(event,maxDispersion,maxZError);
2276  if (!ip) {
2277  ip = FindSimpleIP(event,maxDispersion,maxZError);
2278  if (!ip) {
2279  if (DebugLevel() > 1) AliWarning("No IP for this event found!");
2280  return 0;
2281  }
2282  }
2283  // Good vertex, return it
2284  return ip;
2285 }
2286 //____________________________________________________________________
2288 {
2289  UInt_t evBits = fInputHandler->IsEventSelected();
2290  Bool_t trgOK = (evBits & fOfflineTriggerMask || fOfflineTriggerMask == 0);
2291  if (DebugLevel())
2292  Printf("Trigger bits=0x%08x mask=0x%08x masked=0x%08x -> %s",
2293  evBits, fOfflineTriggerMask, evBits & fOfflineTriggerMask,
2294  trgOK ? "selected" : "rejected");
2295  return trgOK;
2296 }
2297 //____________________________________________________________________
2299 {
2300  static TString centMeth;
2301  if (centMeth.IsNull()) {
2302  centMeth = fCentMethod(3,fCentMethod.Length()-3);
2303  }
2304  AliCentrality* cent = event->GetCentrality();
2305  if (!cent) return -1;
2306 
2307  Double_t centPer = cent->GetCentralityPercentileUnchecked(centMeth);
2308  return centPer;
2309 }
2310 //____________________________________________________________________
2312  Int_t& nTracklets)
2313 {
2314  AliMultSelection* cent =
2315  static_cast<AliMultSelection*>(event->FindListObject("MultSelection"));
2316  if (!cent) {
2317  AliWarning("No centrality in event");
2318  event->GetList()->Print();
2319  return -1;
2320  }
2321  if (fCentIdx < 0) {
2322  TString trkName("SPDTracklets");
2323  for (Int_t i = 0; i < cent->GetNEstimators(); i++) {
2324  AliMultEstimator* e = cent->GetEstimator(i);
2325  if (!e) continue;
2326  if (fCentMethod.EqualTo(e->GetName())) fCentIdx = i;
2327  if (trkName.EqualTo(e->GetName())) fTrkIdx = i;
2328  }
2329  }
2330  if (fCentIdx < 0) {
2331  AliWarningF("Centrality estimator %s not found", fCentMethod.Data());
2332  return -1;
2333  }
2334 
2335  // Use cached index to look up the estimator
2336  AliMultEstimator* estCent = cent->GetEstimator(fCentIdx);
2337  if (!estCent) {
2338  AliWarningF("Centrality estimator %s not available for this event",
2339  fCentMethod.Data());
2340  return -1;
2341  }
2342  fCentEst->Fill(estCent->GetPercentile(), estCent->GetValue());
2343 
2344  // Now look up tracklet value using cached index
2345  if (fTrkIdx < 0) return estCent->GetPercentile();
2346  AliMultEstimator* estTracklets = cent->GetEstimator(fTrkIdx);
2347  if (estTracklets) nTracklets = estTracklets->GetValue();
2348 
2349  return estCent->GetPercentile();
2350 }
2351 
2352 //____________________________________________________________________
2354  Int_t& nTracklets)
2355 {
2356  if (fCentMethod.EqualTo("MB", TString::kIgnoreCase)) {
2357  Printf("MB centrality - not checked");
2358  return 0;
2359  }
2360  Double_t centPer = -1;
2361  if (fCentMethod.BeginsWith("OLD"))
2362  centPer = FindCompatCentrality(event);
2363  else
2364  centPer = FindMultCentrality(event, nTracklets);
2365  const Double_t safety = 1e-3;
2366  if (DebugLevel() > 1) Printf("Read centrality: %f%%", centPer);
2367  if (centPer < -safety) return -2;
2368  if (centPer < +safety) centPer = safety;
2369  else if (centPer > 100-safety) centPer = 100-safety;
2370 
2371  if (fAbsMinCent >= 0 && centPer < fAbsMinCent) {
2372  AliWarningF("Centrality = %f lower than absolute minimum [%f]",
2373  centPer, fAbsMinCent);
2374  return -3;
2375  }
2376  if (centPer < fCentAxis.GetXmin() || centPer > fCentAxis.GetXmax()) {
2377  AliWarningF("Centrality = %f out of range [%f,%f]",
2378  centPer, fCentAxis.GetXmin(), fCentAxis.GetXmax());
2379  return -3;
2380  }
2381  return centPer;
2382 }
2383 
2384 //____________________________________________________________________
2386  Double_t cent,
2387  Double_t ipz)
2388 {
2389  return 1;
2390 }
2391 //____________________________________________________________________
2393  Double_t cent,
2394  Double_t ipz)
2395 {
2396  // We don't check for weights, as we must have them to come this far
2397  // if (!fWeights) {
2398  // AliWarning("No weights defined");
2399  // return 1;
2400  // }
2401  Double_t w = fWeights->LookupWeight(tracklet, cent, ipz, fWeightCorr);
2402  if (tracklet->IsMeasured())
2403  fEtaWeight->Fill(tracklet->GetEta(), cent, w);
2404 
2405  // printf("Looking up weight of tracklet -> %f ", w);
2406  // tracklet->Print();
2407  return w;
2408 }
2409 //____________________________________________________________________
2411 {
2412  Double_t dPhiS = (tracklet->GetDPhi() -
2413  TMath::Sign(fDPhiShift,Double_t(tracklet->GetDPhi())));
2414  UShort_t ret = 0;
2415  ret |= ((tracklet->GetDelta() <= fDeltaCut) ? 0x1 : 0);
2416  ret |= ((dPhiS < fShiftedDPhiCut) ? 0x2 : 0);
2417  // Printf("dPhiS=%f (%f) Delta=%f (%f) -> 0x%x",
2418  // dPhiS, fShiftedDPhiCut,
2419  // tracklet->GetDelta(), fDeltaCut,
2420  // ret);
2421  return ret;
2422 }
2423 
2424 
2425 //____________________________________________________________________
2427  const AliVVertex* ip,
2428  TClonesArray* tracklets)
2429 {
2430  // Figure out which centrality bins to fill
2431  Int_t nAcc = 0;
2432  TIter nextAcc(fCentBins);
2433  CentBin* bin = 0;
2434  TList toRun;
2435  while ((bin = static_cast<CentBin*>(nextAcc()))) {
2436  if (!bin->Accept(cent, ip->GetZ())) continue; // Not in range for this bin
2437  toRun.Add(bin);
2438  nAcc++;
2439  }
2440  // If we have no centrality bins to fill, we return immediately
2441  if (nAcc <= 0) return;
2442 
2443  Double_t ipz = ip->GetZ();
2444  Double_t nBare = 0;
2445  Double_t nMeasured = 0;
2446  Double_t nGenerated = 0;
2447  Double_t nGood = 0;
2448  Double_t nFake = 0;
2449  AliAODTracklet* tracklet = 0;
2450  TIter nextTracklet(tracklets);
2451  while ((tracklet = static_cast<AliAODTracklet*>(nextTracklet()))) {
2452  Double_t weight = LookupWeight(tracklet, cent, ipz);
2453  UShort_t signal = CheckTracklet(tracklet);
2454  if (signal) fEtaPhi->Fill(tracklet->GetEta(), tracklet->GetPhi());
2455  if (tracklet->IsMeasured() && TMath::Abs(tracklet->GetEta()) < 0.7) {
2456  nMeasured += weight;
2457  nBare ++;
2458  if (tracklet->IsCombinatorics()) nFake += weight;
2459  else nGood += weight;
2460  }
2461  else if (tracklet->IsGenerated() &&
2462  !tracklet->IsNeutral() &&
2463  TMath::Abs(tracklet->GetEta()) <= 1) nGenerated ++; // += weight;
2464  TIter nextBin(&toRun);
2465  while ((bin = static_cast<CentBin*>(nextBin()))) {
2466  bin->ProcessTracklet(tracklet, ip->GetZ(), signal, weight);
2467  }
2468  }
2469  fNBareVsFake ->Fill(nBare, nFake);
2470  fNBareVsGood ->Fill(nBare, nGood);
2471  fNTrackletVsFake ->Fill(nMeasured, nFake);
2472  fNTrackletVsGood ->Fill(nMeasured, nGood);
2473  fNGeneratedVsFake->Fill(nGenerated, nFake);
2474  fNGeneratedVsGood->Fill(nGenerated, nGood);
2475 }
2476 
2477 //____________________________________________________________________
2479 {
2480  return fLow >= fHigh;
2481 }
2482 //____________________________________________________________________
2484 {
2485  if (!IsAllBin() && (cent < fLow || cent >= fHigh)) return false;
2486  fCent ->Fill(cent);
2487  fIPz ->Fill(ipz);
2488  fCentIPz->Fill(ipz, cent);
2489  return true;
2490 }
2491 
2492 //____________________________________________________________________
2494  Double_t ipZ,
2495  UShort_t signal,
2496  Double_t weight)
2497 {
2498  TIter next(fSubs);
2499  Histos* h = 0;
2500  if (fDebug > 3) tracklet->Print();
2501  while ((h = static_cast<Histos*>(next())))
2502  h->ProcessTracklet(tracklet, ipZ, signal, weight);
2503 
2504  return true;
2505 }
2506 //____________________________________________________________________
2508  Double_t ipZ,
2509  UShort_t signal,
2510  Double_t weight)
2511 {
2512  if (!fEtaIPz && !fEtaDeltaIPz) return true;
2513  // Get tracklet info
2514  Double_t eta = tracklet->GetEta();
2515  Double_t delta = tracklet->GetDelta();
2516  UChar_t flags = tracklet->GetFlags();
2517  Int_t pdgBin = -1;
2518 
2519  // For debugging
2520  char m[7];
2521  if (fDebug > 3) Bits2String(flags, m);
2522 
2523  // Check the filter mask
2524  if (fMask != 0 && (flags & fMask) == 0) {
2525  if (fDebug > 3) Printf("%14s (0x%02x,----) %6s %7s (0x%02x) ",
2526  GetName(),fMask,"reject",m,flags);
2527  return false;
2528  }
2529 
2530  // If we need the PDG code, get that here
2531  if (fEtaPdgIPz || fEtaPdg || fEtaDeltaPdg)
2532  pdgBin = PdgBin(tracklet->GetParentPdg());
2533 
2534  // If we have the eta-vs-pdg histogram, we should fill before the
2535  // veto, which filters out the neutral particles.
2536  if (fEtaPdg) fEtaPdg->Fill(eta, pdgBin, weight);
2537 
2538  // Check the veto mask
2539  if (fVeto != 0 && (flags & fVeto) != 0) {
2540  if (fDebug > 3) Printf("%14s (----,0x%02x) %6s %7s (0x%02x) ",
2541  GetName(), fVeto, "veto", m, flags);
2542  return false;
2543  }
2544  // Debug message that we accepted tracklet
2545  if (fDebug > 3) Printf("%14s (0x%02x,0x%02x) %6s %7s (0x%02x) ",
2546  GetName(), fMask, fVeto, "accept", m, flags);
2547 
2548  // Fill PDG,eta dependent Delta distributions
2549  if (fEtaDeltaPdg) {
2550  fEtaDeltaPdg->Fill(eta, delta, pdgBin, weight);
2551  if (tracklet->IsCombinatorics())
2552  // Fill both particles
2553  fEtaPdgIPz->Fill(eta,delta,PdgBin(tracklet->GetParentPdg(true)),weight);
2554  }
2555 
2556  // Both reguirements (Delta < cut, dPhi < cut)
2557  if (fEtaIPz && (signal == 0x3)) fEtaIPz->Fill(eta, ipZ, weight);
2558  // Just check dPhi
2559  if (fEtaDeltaIPz && (signal & 0x2)) fEtaDeltaIPz->Fill(eta,delta,ipZ,weight);
2560 
2561  // If we do not need the eta-PDG-IPz or eta-Pt, return now
2562  if (!fEtaPdgIPz && !fEtaPt) return true;
2563 
2564  if (fEtaPdgIPz) {
2565  fEtaPdgIPz->Fill(eta, pdgBin, ipZ, weight);
2566  if (tracklet->IsCombinatorics())
2567  // Fill both particles
2568  fEtaPdgIPz->Fill(eta, PdgBin(tracklet->GetParentPdg(true)), ipZ, weight);
2569  }
2570 
2571  if (fEtaPt) fEtaPt ->Fill(eta, tracklet->GetParentPt(), weight);
2572 
2573  return true;
2574 }
2575 
2576 
2577 //____________________________________________________________________
2578 void
2580 {
2581  Container* results = new Container;
2582  results->SetName(Form("%sResults",GetName()));
2583  results->SetOwner();
2584 
2585  Print("");
2586  fContainer = static_cast<Container*>(GetOutputData(1));
2587  if (!fContainer) {
2588  AliWarning("No sum container found!");
2589  return;
2590  }
2591 
2592  if (!InitCentBins(fContainer)) {
2593  AliWarningF("Failed to initialize centrality bins from %s",
2594  fContainer->GetName());
2595  return;
2596  }
2597 
2598  if (!MasterFinalize(results)) {
2599  AliWarning("Failed to finalize results");
2600  return;
2601  }
2602 
2603  PostData(2, results);
2604 }
2605 
2606 
2607 //____________________________________________________________________
2609 {
2610  fContainer = GetC(parent, fName);
2611  fCent = GetH1(fContainer, "cent");
2612  fIPz = GetH1(fContainer, "ipz");
2613  fCentIPz = GetP1(fContainer, "centIpz");
2614  if (!fContainer || !fCent || !fIPz) return false;
2615  TIter next(fSubs);
2616  Histos* h = 0;
2617  while ((h = static_cast<Histos*>(next())))
2618  if (!h->FinalizeInit(fContainer)) return false;
2619  return true;
2620 }
2621 //____________________________________________________________________
2623 {
2624  fContainer = GetC(parent, fName);
2625  fEtaIPz = GetH2(fContainer, "etaIPz", false); // No complaints
2626  fEtaDeltaIPz = GetH3(fContainer, "etaDeltaIPz", false);
2627  fEtaDeltaPdg = GetH3(fContainer, "etaDeltaPdg", false);
2628  fEtaPdgIPz = GetH3(fContainer, "etaPdgIPz", false);
2629  fEtaPdg = GetH2(fContainer, "etaPdg", false);
2630  fEtaPt = GetH2(fContainer, "etaPt", false);
2631  if (GetName()[0] == 'm' && GetC(parent,"generated")) {
2632  // Fix up titles
2633  if (fEtaIPz)
2634  fEtaIPz->SetTitle(Form("%s'", fEtaIPz->GetTitle()));
2635  if (fEtaDeltaIPz)
2636  fEtaDeltaIPz->SetTitle(Form("#Delta_{%s}", fEtaIPz->GetTitle()));
2637  }
2638  return (fContainer != 0); // && fEtaIPz != 0);
2639 }
2640 
2641 //____________________________________________________________________
2643 {
2644  // Make copies of histograms and store
2645  fIPz = static_cast<TH1*>(CloneAndAdd(results,GetH1(fContainer,"ipz")));
2646  fCent = static_cast<TH1*>(CloneAndAdd(results,GetH1(fContainer,"cent")));
2647  fStatus = static_cast<TH1*>(CloneAndAdd(results,GetH1(fContainer,"status")));
2648  fEtaPhi = static_cast<TH2*>(CloneAndAdd(results,GetH2(fContainer,"etaPhi")));
2649  fCentTracklets =
2650  static_cast<TProfile*>(CloneAndAdd(results,GetP(fContainer,
2651  "centTracklets")));
2652  fCentEst =
2653  static_cast<TProfile*>(CloneAndAdd(results,GetP(fContainer,
2654  "centEstimator")));
2655  Double_t nEvents = fIPz->GetEntries();
2656  Printf("Event summary:");
2657  for (Int_t i = 1; i <= fStatus->GetNbinsX(); i++)
2658  Printf(" %10d %s",
2659  Int_t(fStatus->GetBinContent(i)),
2660  fStatus->GetXaxis()->GetBinLabel(i));
2661  for (Int_t i = 1; i <= fCent->GetNbinsX(); i++)
2662  Printf(" %6.2f-%6.2f%%: %d",
2663  fCent->GetXaxis()->GetBinLowEdge(i),
2664  fCent->GetXaxis()->GetBinUpEdge(i),
2665  Int_t(fCent->GetBinContent(i)));
2666 
2667 
2668  fIPz ->Scale(1./nEvents);
2669  fCent ->Scale(1./fCent->GetEntries());
2670  fStatus->Scale(1./fStatus->GetBinContent(1));
2671  fEtaPhi->Scale(1./nEvents);
2672 
2673  if (fTailMax < 0) fTailMax = fMaxDelta;
2674  TIter next(fCentBins);
2675  CentBin* bin = 0;
2676  while ((bin = static_cast<CentBin*>(next()))) {
2677  if (!bin->MasterFinalize(results, 0, fTailDelta, fTailMax)) {
2678  AliWarningF("Failed to finalize %s", bin->GetName());
2679  return false;
2680  }
2681  }
2682  return true;
2683 }
2684 //____________________________________________________________________
2686 {
2687  if (!AliTrackletAODMCdNdeta::MasterFinalize(results)) return false;
2688 
2689  TObject* o = fContainer->FindObject("etaWeight");
2690  if (o && o->IsA()->InheritsFrom(TH1::Class())) {
2691  TH1* h = static_cast<TH1*>(o->Clone());
2692  h->SetDirectory(0);
2693  results->Add(h);
2694  }
2695  else {
2696  AliWarningF("Object %p (etaWeight) is not a TH1 or not found",o);
2697  }
2698  o = fContainer->FindObject("weightCorr");
2699  if (o && o->IsA()->InheritsFrom(TH1::Class())) {
2700  TH1* h = static_cast<TH1*>(o->Clone());
2701  h->SetDirectory(0);
2702  results->Add(h);
2703  }
2704  else {
2705  AliWarningF("Object %p (weightCorr) is not a TH1 or not found",o);
2706  }
2707  if (!fWeights) return true;
2708 
2709  if (!fWeights->Retrieve(fContainer)) return false;
2710  fWeights->Store(results);
2711 
2712  return true;
2713 }
2714 
2715 
2716 //____________________________________________________________________
2718  TH1* ,
2719  Double_t tailDelta,
2720  Double_t tailMax)
2721 {
2722  Container* result = new Container;
2723  result->SetName(fName);
2724  result->SetOwner(true);
2725  parent->Add(result);
2726 
2727  Double_t nEvents = fIPz->GetEntries();
2728  // Copy ipZ histogram and scale by number of events
2729  TH1* ipZ = static_cast<TH1*>(CloneAndAdd(result, fIPz));
2730  ipZ->Scale(1./nEvents);
2731 
2732  TH1* cent = static_cast<TH1*>(CloneAndAdd(result, fCent));
2733  cent->Scale(1./nEvents);
2734 
2735  CloneAndAdd(result, fCentIPz);
2736 
2737  Container* measCont = 0;
2738  Container* genCont = 0;
2739  TIter next(fSubs);
2740  Histos* h = 0;
2741  while ((h = static_cast<Histos*>(next()))) {
2742  if (h->GetMask() != AliAODTracklet::kInjection &&
2744  // For the anything but injection or MC-labels, we just finalize
2745  if (!h->MasterFinalize(result, fIPz, tailDelta, tailMax)) {
2746  AliWarningF("Failed to finalize %s/%s", GetName(), h->GetName());
2747  return false;
2748  }
2749  if (h->GetMask() == AliAODTracklet::kGenerated)
2750  // Set MC truth container
2751  genCont = GetC(result, h->GetName());
2752  if (h == fMeasured)
2753  measCont = GetC(result, h->GetName());
2754  continue;
2755  }
2756  if (!EstimateBackground(result, measCont, genCont, h, tailDelta, tailMax)) {
2757  AliWarningF("Failed to estimate Bg in %s/%s", GetName(), h->GetName());
2758  return false;
2759  }
2760  }
2761  return true;
2762 }
2763 
2764 //____________________________________________________________________
2765 Bool_t
2767  Int_t nEvents)
2768 {
2769  // Scale distribution of PDGs to number of events and bin width
2770  if (!fEtaPdg) return true;
2771 
2772 
2773  TH2* etaPdg = static_cast<TH2*>(fEtaPdg->Clone());
2774  etaPdg->SetDirectory(0);
2775  etaPdg->Scale(1. / nEvents, "width");
2776  result->Add(etaPdg);
2777 
2778  // Loop over PDG types and create 2D and 1D distributions
2779  TAxis* yaxis = etaPdg->GetYaxis();
2780  Int_t first = yaxis->GetFirst();
2781  Int_t last = yaxis->GetLast();
2782  THStack* pdgs = new THStack("all","");
2783  THStack* toPion = new THStack("toPion", "");
2784  THStack* toAll = new THStack("toAll", "");
2785  TH1* pion = 0;
2786  Container* pdgOut = new Container();
2787  pdgOut->SetName("mix");
2788  result->Add(pdgOut);
2789  pdgOut->Add(pdgs);
2790  pdgOut->Add(toPion);
2791  pdgOut->Add(toAll);
2792 
2793  TH1* all = static_cast<TH1*>(etaPdg->ProjectionX("total",
2794  1,etaPdg->GetNbinsY()));
2795  all->SetDirectory(0);
2796  all->SetName("total");
2797  all->SetTitle("All");
2798  all->SetYTitle(Form("d#it{N}_{%s}/d#eta", "all"));
2799  all->SetFillColor(kBlack);
2800  all->SetMarkerColor(kBlack);
2801  all->SetLineColor(kBlack);
2802  all->SetMarkerStyle(20);
2803  all->SetFillStyle(0);
2804  all->SetFillColor(0);
2805  all->Reset();
2806  pdgs->Add(all);
2807  for (Int_t i = 1; i <= etaPdg->GetNbinsY(); i++) {
2808  Int_t pdg = TString(yaxis->GetBinLabel(i)).Atoi();
2809  TString nme;
2810  Style_t sty;
2811  Color_t col;
2812  PdgAttr(pdg, nme, col, sty);
2813  if (pdg < 0) pdg = 0;
2814  if (pdg == 22) continue; // Ignore photons
2815 
2816  TH1* h1 = static_cast<TH1*>(etaPdg->ProjectionX(Form("h%d", pdg),i,i));
2817  if (h1->GetEntries() <= 0) continue; // Do not store if empty
2818  h1->SetDirectory(0);
2819  h1->SetName(Form("eta_%d", pdg));
2820  h1->SetTitle(nme);
2821  h1->SetYTitle(Form("d#it{N}_{%s}/d#eta", nme.Data()));
2822  h1->SetFillColor(col);
2823  h1->SetMarkerColor(col);
2824  h1->SetLineColor(col);
2825  h1->SetMarkerStyle(sty);
2826  h1->SetFillStyle(0);
2827  h1->SetFillColor(0);
2828  h1->SetBinContent(0,0);
2829  all->Add(h1);
2830  pdgs->Add(h1);
2831  switch (pdg) {
2832  case 321: h1->SetBinContent(0, 0.15); break; // PRC88,044910
2833  case 2212: h1->SetBinContent(0, 0.05); break; // PRC88,044910
2834  case 310: h1->SetBinContent(0, 0.075); break; // PRL111,222301
2835  case 3122: h1->SetBinContent(0, 0.018); break; // PRL111,222301
2836  case 3212: h1->SetBinContent(0, 0.0055); break; // NPA904,539
2837  case 3322: h1->SetBinContent(0, 0.005); break; // PLB734,409
2838  case 211: h1->SetBinContent(0, 1); break; // it self
2839  default: h1->SetBinContent(0, -1); break; // Unknown
2840  }
2841  pdgOut->Add(h1);
2842 
2843  if (pdg == 211) pion = h1;
2844  }
2845  if (!pdgs->GetHists()) return true;
2846 
2847  TIter next(pdgs->GetHists());
2848  TH1* tmp = 0;
2849  Double_t rmin = +1e9;
2850  Double_t rmax = -1e9;
2851  while ((tmp = static_cast<TH1*>(next()))) {
2852  if (tmp == all) continue;
2853  // Calculate ratio to all
2854  TH1* rat = static_cast<TH1*>(tmp->Clone());
2855  rat->Divide(all);
2856  rat->SetDirectory(0);
2857  rat->SetTitle(Form("%s / all", tmp->GetTitle()));
2858  toAll->Add(rat);
2859 
2860  if (tmp == pion) continue;
2861  Double_t r276 = tmp->GetBinContent(0);
2862  if (r276 < 0 || r276 >= 1) continue;
2863 
2864  // Calulate ratio to pions
2865  rat = static_cast<TH1*>(tmp->Clone());
2866  rat->Divide(pion);
2867  rat->SetTitle(Form("%s / %s", tmp->GetTitle(), pion->GetTitle()));
2868  rat->SetDirectory(0);
2869 
2870  TGraphErrors* g = new TGraphErrors(1);
2871  g->SetName(Form("%s_2760", rat->GetName()));
2872  g->SetTitle(Form("%s in #sqrt{s_{NN}}=2.76TeV", rat->GetTitle()));
2873  g->SetPoint(0,0,r276);
2874  g->SetPointError(0,.5,0);
2875  g->SetLineColor(rat->GetLineColor());
2876  g->SetLineStyle(rat->GetLineStyle());
2877  g->SetMarkerColor(rat->GetMarkerColor());
2878  g->SetMarkerStyle(rat->GetMarkerStyle());
2879  g->SetMarkerSize(1.5*rat->GetMarkerSize());
2880  rat->GetListOfFunctions()->Add(g,"p");
2881  rat->SetMaximum(TMath::Max(rat->GetMaximum(),r276));
2882  rat->SetMinimum(TMath::Min(rat->GetMinimum(),r276));
2883  rmin = TMath::Min(rat->GetMinimum(),rmin);
2884  rmax = TMath::Max(rat->GetMaximum(),rmax);
2885 
2886  toPion->Add(rat);
2887  }
2888  // toPion->SetMinimum(rmin);
2889  toPion->SetMaximum(1.1*rmax);
2890 
2891  return true;
2892 }
2893 
2894 //____________________________________________________________________
2895 Bool_t
2897  Int_t nEvents,
2898  Double_t tailDelta,
2899  Double_t tailMax,
2900  const TString& pre,
2901  const TString& tit)
2902 {
2903  Container* pdgOut = new Container();
2904  pdgOut->SetName(pre);
2905  result->Add(pdgOut);
2906  TAxis* zaxis = fEtaDeltaPdg->GetZaxis();
2907  Int_t first = zaxis->GetFirst();
2908  Int_t last = zaxis->GetLast();
2909  Bool_t isMid = TString(pre).EqualTo("mid");
2910  Int_t nX = fEtaDeltaPdg->GetNbinsX();
2911  Int_t nY = fEtaDeltaPdg->GetNbinsY();
2912  Int_t nZ = fEtaDeltaPdg->GetNbinsZ();
2913  Int_t l1 = fEtaDeltaPdg->GetXaxis()->FindBin(-1);
2914  Int_t r1 = fEtaDeltaPdg->GetXaxis()->FindBin(+1);
2915  Int_t b1 = (isMid ? l1 : 1);
2916  Int_t b2 = (isMid ? r1 : l1-1);
2917  Int_t b3 = (isMid ? -1 : r1+1);
2918  Int_t b4 = (isMid ? -1 : nX);
2919  THStack* stack = new THStack("all", tit);
2920  pdgOut->Add(stack);
2921 
2922  TH1* total = fEtaDeltaPdg->ProjectionY("totalMid",1, 1,1, nZ);
2923  total->SetDirectory(0);
2924  total->SetName("total");
2925  total->SetTitle("Total");
2926  total->SetYTitle(Form("d#it{N}_{%s}/d#Delta", "total"));
2927  total->SetFillColor(kBlack);
2928  total->SetMarkerColor(kBlack);
2929  total->SetLineColor(kBlack);
2930  total->SetMarkerStyle(24);
2931  total->SetFillStyle(0);
2932  total->SetFillColor(0);
2933  total->Reset();
2934  stack->Add(total);
2935 
2936  THStack* ratios = new THStack("ratios","");
2937  TH1* ratioSig = Make1D(0, "ratioSig", "Ratio to all", kGreen+1, 24, *zaxis);
2938  TH1* ratioBg = Make1D(0, "ratioBg", "Ratio to all", kRed+1, 25, *zaxis);
2939  ratioSig->GetXaxis()->LabelsOption("v");
2940  ratioBg ->GetXaxis()->LabelsOption("v");
2941  ratioSig->SetFillColor(kGreen+1);
2942  ratioBg ->SetFillColor(kRed +1);
2943  ratioSig->SetFillStyle(1001);
2944  ratioBg ->SetFillStyle(1001);
2945  ratioSig->SetBarWidth(0.4);
2946  ratioBg ->SetBarWidth(0.4);
2947  ratioSig->SetBarOffset(0.05);
2948  ratioBg ->SetBarOffset(0.55);
2949  ratios->Add(ratioSig, "bar0 e text30");
2950  ratios->Add(ratioBg, "bar0 e text30");
2951  pdgOut->Add(ratios);
2952 
2953  THStack* rows = new THStack("rows","");
2954  TH1* rowSig = new TH1D("rowSig","row",6,0,6);
2955  rowSig->SetDirectory(0);
2956  rowSig->SetYTitle("Fraction of signal");
2957  rowSig->GetXaxis()->SetBinLabel(1, "K^{0}_{S}");
2958  rowSig->GetXaxis()->SetBinLabel(2, "K^{#pm}");
2959  rowSig->GetXaxis()->SetBinLabel(3, "#Lambda");
2960  rowSig->GetXaxis()->SetBinLabel(4, "#Xi");
2961  rowSig->GetXaxis()->SetBinLabel(5, "#Sigma");
2962  rowSig->GetXaxis()->SetBinLabel(6, "Other");
2963  rowSig->GetXaxis()->LabelsOption("v");
2964  rowSig->SetMaximum(100);
2965  rowSig->SetMinimum(0);
2966  rowSig->SetLineColor(kGreen+1);
2967  rowSig->SetMarkerColor(kGreen+1);
2968  rowSig->SetMarkerStyle(24);
2969  TH1* rowBg = static_cast<TH1*>(rowSig->Clone("rowBg"));
2970  rowBg->SetDirectory(0);
2971  rowBg->SetYTitle("Fraction of background");
2972  rowBg->SetLineColor(kRed+1);
2973  rowBg->SetMarkerColor(kRed+1);
2974  rowBg->SetMarkerStyle(25);
2975  rowSig->SetFillStyle(1001);
2976  rowBg ->SetFillColor(1001);
2977  rowSig->SetFillColor(kGreen+1);
2978  rowBg ->SetFillColor(kRed +1);
2979  rowSig->SetBarWidth(0.4);
2980  rowBg ->SetBarWidth(0.4);
2981  rowSig->SetBarOffset(0.05);
2982  rowBg ->SetBarOffset(0.55);
2983  rows->Add(rowSig, "bar0 e text30");
2984  rows->Add(rowBg, "bar0 e text30");
2985  pdgOut->Add(rows);
2986 
2987  for (Int_t i = 1; i <= zaxis->GetNbins(); i++) {
2988  Int_t pdg = TString(zaxis->GetBinLabel(i)).Atoi();
2989  TString nme;
2990  Style_t sty;
2991  Color_t col;
2992  Int_t ipdg = pdg;
2993  PdgAttr(pdg, nme, col, sty);
2994  if (pdg < 0) pdg = 0;
2995  ratioSig->GetXaxis()->SetBinLabel(i, nme);
2996  ratioBg ->GetXaxis()->SetBinLabel(i, nme);
2997  if (pdg == 22) continue; // Ignore photons
2998 
2999 
3000  TH1* h1 = fEtaDeltaPdg->ProjectionY(Form("%d", pdg), b1, b2, i,i);
3001  h1->SetDirectory(0);
3002  if (b3 < b4) {
3003  TH1* h2 = fEtaDeltaPdg->ProjectionY(Form("t%d", pdg), b3, b4, i,i);
3004  h2->SetDirectory(0);
3005  h1->Add(h2);
3006  delete h2;
3007  }
3008  h1->SetUniqueID(i); // Store bin number
3009  if (h1->GetEntries() <= 0) continue; // Do not store if empty
3010  h1->SetName(Form("%d", pdg));
3011  h1->SetTitle(nme);
3012  h1->SetYTitle(Form("d#it{N}_{%s}/d#Delta", nme.Data()));
3013  h1->SetFillColor(col);
3014  h1->SetMarkerColor(col);
3015  h1->SetLineColor(col);
3016  h1->SetMarkerStyle(sty);
3017  h1->SetFillStyle(0);
3018  h1->SetFillColor(0);
3019  h1->SetBinContent(0,ipdg);
3020  total->Add(h1);
3021  stack->Add(h1);
3022  }
3023  total->SetBinContent(0,0);
3024  total->Scale(1./nEvents);
3025  Double_t deltaCut = 1.5;
3026  Double_t eTot, iTot = Integrate(total, 0, tailMax, eTot);
3027  Double_t eSig, iSig = Integrate(total, 0, deltaCut, eSig);
3028  Double_t eBg, iBg = Integrate(total, tailDelta, tailMax, eBg);
3029 
3030  TH1* totInt = new TH1D("totalIntegrals", "Integrals of total:", 3, 0, 3);
3031  totInt->SetDirectory(0);
3032  totInt->GetXaxis()->SetBinLabel(1,Form("0#minus%4.1f", tailMax));
3033  totInt->GetXaxis()->SetBinLabel(2,Form("0#minus%3.1f",deltaCut));
3034  totInt->GetXaxis()->SetBinLabel(3,Form("%3.1f#minus%4.1f",tailDelta,tailMax));
3035  totInt->SetBinContent(1,iTot); totInt->SetBinError(1, eTot);
3036  totInt->SetBinContent(2,iSig); totInt->SetBinError(2, eSig);
3037  totInt->SetBinContent(3,iBg); totInt->SetBinError(3, eBg);
3038  pdgOut->Add(totInt);
3039 
3040  if (!stack->GetHists()) return true;
3041 
3042  Double_t sigOth = 0, eSigOth = 0, bgOth = 0, eBgOth = 0;
3043  Double_t sigK0s = 0, eSigK0s = 0, bgK0s = 0, eBgK0s = 0;
3044  Double_t sigKpm = 0, eSigKpm = 0, bgKpm = 0, eBgKpm = 0;
3045  Double_t sigLam = 0, eSigLam = 0, bgLam = 0, eBgLam = 0;
3046  Double_t sigXi0 = 0, eSigXi0 = 0, bgXi0 = 0, eBgXi0 = 0;
3047  Double_t sigSig = 0, eSigSig = 0, bgSig = 0, eBgSig = 0;
3048  TIter next(stack->GetHists());
3049  TH1* tmp = 0;
3050  while ((tmp = static_cast<TH1*>(next()))) {
3051  // Calculate ratio to all
3052  // Scale(tmp, iTot, eTot);
3053  if (tmp == total) continue;
3054  Int_t pdg = tmp->GetBinContent(0); tmp->SetBinContent(0,0);
3055  tmp->Scale(1./ nEvents);
3056  Double_t eiSig, iiSig = Integrate(tmp, 0, deltaCut, eiSig);
3057  Double_t eiBg, iiBg = Integrate(tmp, tailDelta, tailMax, eiBg);
3058  Double_t erS, rS = RatioE(iiSig, eiSig, iSig, eSig, erS);
3059  Double_t erB, rB = RatioE(iiBg, eiBg, iBg, eBg, erB);
3060  Int_t bin = tmp->GetUniqueID();
3061 #if 0
3062  Printf(" %10s(%4d,%3d) "
3063  "signal=%6.1f/%6.1f=%5.3f+/-%5.3f bg=%6.1f/%6.1f=%5.3f+/-%5.3f",
3064  tmp->GetTitle(), pdg, bin, iiSig, iSig, rS, erS, iiBg, iBg, rB, erB);
3065 #endif
3066  ratioSig->SetBinContent(bin, rS);
3067  ratioSig->SetBinError (bin, erS);
3068  ratioBg ->SetBinContent(bin, rB);
3069  ratioBg ->SetBinError (bin, erB);
3070 
3071  switch (pdg) {
3072  case 321: // K^{+}
3073  sigKpm += rS; eSigKpm += erS*erS; bgKpm += rB; eBgKpm += erB*erB; break;
3074  case 310: // K^{0}_{S}
3075  sigK0s += rS; eSigK0s += erS*erS; bgK0s += rB; eBgK0s += erB*erB; break;
3076  case 3112: // #Sigma^{-}
3077  case 3222: // #Sigma^{+}
3078  // case 3114: // #Sigma^{*-}
3079  // case 3224: // #Sigma^{*+}
3080  // case 4214: // #Sigma^{*+}_{c}
3081  // case 4224: // #Sigma^{*++}_{c}
3082  // case 3212: // #Sigma^{0}
3083  // case 4114: // #Sigma^{*0}_{c}
3084  // case 3214: // #Sigma^{*0}
3085  sigSig += rS; eSigSig += erS*erS; bgSig += rB; eBgSig += erB*erB; break;
3086  case 3312: // #Xi^{-}
3087  // case 3314: // #Xi^{*-}
3088  // case 3324: // #Xi^{*0}
3089  // case 4132: // #Xi^{0}_{c}
3090  // case 4314: // #Xi^{*0}_{c}
3091  sigXi0 += rS; eSigXi0 += erS*erS; bgXi0 += rB; eBgXi0 += erB*erB; break;
3092  // case 4122: // #Lambda^{+}_{c}
3093  case 3122: // #Lambda
3094  sigLam += rS; eSigLam += erS*erS; bgLam += rB; eBgLam += erB*erB; break;
3095  default:
3096  sigOth += rS; eSigOth += erS*erS; bgOth += rB; eBgOth += erB*erB; break;
3097  }
3098  } // while(tmp)
3099  rowSig->SetBinContent(1,100*sigK0s);
3100  rowSig->SetBinError (1,100*TMath::Sqrt(eSigK0s));
3101  rowSig->SetBinContent(2,100*sigKpm);
3102  rowSig->SetBinError (2,100*TMath::Sqrt(eSigKpm));
3103  rowSig->SetBinContent(3,100*sigLam);
3104  rowSig->SetBinError (3,100*TMath::Sqrt(eSigLam));
3105  rowSig->SetBinContent(4,100*sigXi0);
3106  rowSig->SetBinError (4,100*TMath::Sqrt(eSigXi0));
3107  rowSig->SetBinContent(5,100*sigSig);
3108  rowSig->SetBinError (5,100*TMath::Sqrt(eSigSig));
3109  rowSig->SetBinContent(6,100*sigOth);
3110  rowSig->SetBinError (6,100*TMath::Sqrt(eSigOth));
3111  rowBg ->SetBinContent(1,100* bgK0s);
3112  rowBg ->SetBinError (1,100*TMath::Sqrt( eBgK0s));
3113  rowBg ->SetBinContent(2,100* bgKpm);
3114  rowBg ->SetBinError (2,100*TMath::Sqrt( eBgKpm));
3115  rowBg ->SetBinContent(3,100* bgLam);
3116  rowBg ->SetBinError (3,100*TMath::Sqrt( eBgLam));
3117  rowBg ->SetBinContent(4,100* bgXi0);
3118  rowBg ->SetBinError (4,100*TMath::Sqrt( eBgXi0));
3119  rowBg ->SetBinContent(5,100* bgSig);
3120  rowBg ->SetBinError (5,100*TMath::Sqrt( eBgSig));
3121  rowBg ->SetBinContent(6,100* bgOth);
3122  rowBg ->SetBinError (6,100*TMath::Sqrt( eBgOth));
3123 
3124  return true;
3125 }
3126 
3127 //____________________________________________________________________
3128 Bool_t
3130  Int_t nEvents,
3131  Double_t tailDelta,
3132  Double_t tailMax)
3133 {
3134  // Scale distribution of PDGs to number of events and bin width
3135  if (!fEtaDeltaPdg) return true;
3136 
3137  Container* pdgOut = new Container();
3138  pdgOut->SetName("specie");
3139  result->Add(pdgOut);
3140 
3141  ProjectEtaDeltaPdgPart(pdgOut,
3142  nEvents,
3143  tailDelta,
3144  tailMax,
3145  "mid",
3146  "|#eta|<1");
3147  ProjectEtaDeltaPdgPart(pdgOut,
3148  nEvents,
3149  tailDelta,
3150  tailMax,
3151  "fwd",
3152  "|#eta|>1");
3153 }
3154 
3155 //____________________________________________________________________
3156 Bool_t
3158  TH1* ipz,
3159  const TString& shn)
3160 {
3161  // If we have the PDG distributions, we project on eta,IPz, and then
3162  // average over IPz to get dNpdg/deta.
3163  if (!fEtaPdgIPz) return true;
3164 
3165  // Scale each vertex range by number of events in that range
3166  TH3* etaPdgIPz = ScaleToIPz(fEtaPdgIPz, ipz, false);
3167  result->Add(etaPdgIPz);
3168 
3169  // Loop over PDG types and create 2D and 1D distributions
3170  TAxis* yaxis = etaPdgIPz->GetYaxis();
3171  Int_t first = yaxis->GetFirst();
3172  Int_t last = yaxis->GetLast();
3173  THStack* pdgs = new THStack("all","");
3174  THStack* ratios = new THStack("toPion", "");
3175  TH1* pion = 0;
3176  Container* pdgOut = new Container();
3177  pdgOut->SetName("types");
3178  result->Add(pdgOut);
3179  pdgOut->Add(pdgs);
3180  pdgOut->Add(ratios);
3181  for (Int_t i = 1; i <= etaPdgIPz->GetNbinsY(); i++) {
3182  yaxis->SetRange(i,i);
3183 
3184  Int_t pdg = TString(yaxis->GetBinLabel(i)).Atoi();
3185  TString nme;
3186  Style_t sty;
3187  Color_t col;
3188  PdgAttr(pdg, nme, col, sty);
3189  if (pdg < 0) pdg = 0;
3190 
3191  TH2* h2 = static_cast<TH2*>(etaPdgIPz->Project3D("zx e"));
3192  if (h2->GetEntries() <= 0) continue; // Do not store if empty
3193  h2->SetDirectory(0);
3194  h2->SetName(Form("etaIPz_%d", pdg));
3195  h2->SetTitle(Form("%s#rightarrowX_{%s}", nme.Data(), shn.Data()));
3196  h2->SetYTitle(Form("d#it{N}^{2}_{%s#rightarrow%s}/"
3197  "(d#etadIP_{#it{z}})",
3198  nme.Data(), shn.Data()));
3199  h2->SetFillColor(col);
3200  h2->SetMarkerColor(col);
3201  h2->SetLineColor(col);
3202  pdgOut->Add(h2);
3203 
3204  TH1* h1 = AverageOverIPz(h2, Form("eta_%d",pdg), 1, ipz, 0, false);
3205  h1->SetDirectory(0);
3206  h1->SetYTitle(Form("d#it{N}_{%s#rightarrow%s}/d#eta",
3207  nme.Data(), shn.Data()));
3208  h1->SetMarkerStyle(sty);
3209  pdgs->Add(h1);
3210 
3211  if (pdg == 211) pion = h1;
3212  }
3213  if (!pdgs->GetHists()) {
3214  yaxis->SetRange(first, last);
3215  return true;
3216  }
3217 
3218  TIter next(pdgs->GetHists());
3219  TH1* tmp = 0;
3220  while ((tmp = static_cast<TH1*>(next()))) {
3221  if (tmp == pion) continue;
3222  TH1* rat = static_cast<TH1*>(tmp->Clone());
3223  rat->Divide(pion);
3224  rat->SetDirectory(0);
3225  ratios->Add(rat);
3226  }
3227 
3228  yaxis->SetRange(first, last);
3229 }
3230 
3231 //____________________________________________________________________
3232 Bool_t
3234  TH1* ipz,
3235  Double_t tailDelta,
3236  Double_t tailMax)
3237 {
3238  Container* result = new Container;
3239  result->SetName(fName);
3240  result->SetOwner(true);
3241  parent->Add(result);
3242 
3243  // Get the number of events
3244  Double_t nEvents = ipz->GetEntries();
3245 
3246  // Scale each vertex range by number of events in that range
3247  TH2* etaIPz = 0;
3248  if (fEtaIPz) {
3249  etaIPz = ScaleToIPz(fEtaIPz, ipz);
3250  result->Add(etaIPz);
3251  }
3252  // Scale distribution of Pt to number of events and bin width
3253  if (fEtaPt) {
3254  TH2* etaPt = static_cast<TH2*>(fEtaPt->Clone());
3255  etaPt->SetDirectory(0);
3256  etaPt->Scale(1. / nEvents, "width");
3257  result->Add(etaPt);
3258  }
3259  // Short-hand-name
3260  TString shn(etaIPz ? etaIPz->GetTitle() : "X");
3261 
3262  ProjectEtaPdg (result, nEvents);
3263  ProjectEtaDeltaPdg(result, nEvents, tailDelta, tailMax);
3264  ProjectEtaPdgIPz (result, ipz, shn);
3265 
3266  // If we do not have eta vs Delta, just return
3267  if (!fEtaDeltaIPz) return true;
3268 
3269  // Normalize delta distribution to integral number of events
3270  // static_cast<TH3*>(CloneAndAdd(result, fEtaDeltaIPz));
3271  TH3* etaDeltaIPz = ScaleToIPz(fEtaDeltaIPz, ipz, false);
3272  // ipz->GetEntries()>1000);
3273  result->Add(etaDeltaIPz);
3274 
3275  // Make 2D projection to eta,Delta
3276  TH2* etaDelta = ProjectEtaDelta(etaDeltaIPz);
3277  result->Add(etaDelta);
3278 
3279  // Make projection of delta
3280  TH1* delta = ProjectDelta(etaDelta);
3281  result->Add(delta);
3282 
3283  // PArameters of integrals
3284  Double_t maxDelta = etaDeltaIPz->GetYaxis()->GetXmax();
3285  Int_t lowBin = etaDeltaIPz->GetYaxis()->FindBin(tailDelta);
3286  Int_t sigBin = etaDeltaIPz->GetYaxis()->FindBin(1.5);
3287  Int_t highBin = TMath::Min(etaDeltaIPz->GetYaxis()->FindBin(tailMax),
3288  etaDeltaIPz->GetYaxis()->GetNbins());
3289 
3290 
3291  TH1* etaDeltaTail = etaDelta->ProjectionX("etaDeltaTail");
3292  etaDeltaTail ->SetDirectory(0);
3293  etaDeltaTail ->Reset();
3294  etaDeltaTail ->SetTitle(Form("#scale[.7]{#int}_{%3.1f}^{%4.1f}"
3295  "d%s d#it{N}/d%s",
3296  tailDelta,maxDelta,
3297  etaDeltaIPz->GetTitle(),
3298  etaDeltaIPz->GetTitle()));
3299  etaDeltaTail ->SetYTitle("#scale[.7]{#int}_{tail}d#Delta d#it{N}/d#Delta");
3300 
3301  TH2* etaIPzDeltaTail = static_cast<TH2*>(etaDeltaIPz->Project3D("zx e"));
3302  etaIPzDeltaTail->SetName("etaIPzDeltaTail");
3303  etaIPzDeltaTail->SetDirectory(0);
3304  etaIPzDeltaTail->Reset();
3305  etaIPzDeltaTail->SetTitle(etaDeltaTail->GetTitle());
3306  etaIPzDeltaTail->SetZTitle(etaDelta->GetYaxis()->GetTitle());
3307  TH2* etaIPzDeltaMain = static_cast<TH2*>(etaDeltaIPz->Project3D("zx e"));
3308  etaIPzDeltaMain->SetName("etaIPzDeltaMain");
3309  etaIPzDeltaMain->SetDirectory(0);
3310  etaIPzDeltaMain->Reset();
3311  etaIPzDeltaMain->SetTitle(etaDeltaTail->GetTitle());
3312  etaIPzDeltaMain->SetZTitle(etaDelta->GetYaxis()->GetTitle());
3313  // Loop over eta
3314  Double_t intg = 0, eintg = 0;
3315  for (Int_t i = 1; i <= etaDeltaTail->GetNbinsX(); i++) {
3316  // Integrate over Delta
3317  intg = etaDelta->IntegralAndError(i, i, lowBin, highBin, eintg);
3318  etaDeltaTail->SetBinContent(i, intg);
3319  etaDeltaTail->SetBinError (i, eintg);
3320  // Loop over IPz
3321  for (Int_t j = 1; j <= etaIPzDeltaTail->GetNbinsY(); j++) {
3322  // Integrate over Delta
3323  intg = etaDeltaIPz->IntegralAndError(i,i,lowBin,highBin,j,j,eintg);
3324  etaIPzDeltaTail->SetBinContent(i, j, intg);
3325  etaIPzDeltaTail->SetBinError (i, j, eintg);
3326 
3327  intg = etaDeltaIPz->IntegralAndError(i,i,1,sigBin,j,j,eintg);
3328  etaIPzDeltaMain->SetBinContent(i, j, intg);
3329  etaIPzDeltaMain->SetBinError (i, j, eintg);
3330  }
3331  }
3332  result->Add(etaIPzDeltaTail);
3333  result->Add(etaIPzDeltaMain);
3334  result->Add(etaDeltaTail);
3335 
3336  // Integrate full tail
3337  intg = etaDeltaIPz->IntegralAndError(1,etaDeltaIPz->GetNbinsX(),
3338  lowBin, highBin,
3339  1,etaDeltaIPz->GetNbinsZ(),
3340  eintg);
3341  result->Add(new TParameter<double>("deltaTail", intg));
3342  result->Add(new TParameter<double>("deltaTailError", eintg));
3343 
3344  // Some consistency checks:
3345  if (fDebug > 1) {
3346  Printf("%10s: Integral over eta,IPz: %9.4f +/- %9.4f",
3347  GetName(), intg, eintg);
3348  intg = etaDelta->IntegralAndError(1,etaDeltaIPz->GetNbinsX(),
3349  lowBin, highBin,
3350  eintg);
3351  Printf("%10s: Integral over eta: %9.4f +/- %9.4f",
3352  GetName(), intg, eintg);
3353  intg = delta->IntegralAndError(lowBin, highBin, eintg);
3354  Printf("%10s: Integral: %9.4f +/- %9.4f",
3355  GetName(), intg, eintg);
3356  }
3357 
3358  return true;
3359 }
3360 
3361 //____________________________________________________________________
3363  Container* measCont,
3364  Container* genCont,
3365  Histos* h,
3366  Double_t tailCut,
3367  Double_t tailMax)
3368 {
3369  if (!h || !measCont) {
3370  AliWarningF("No sub-histos or measured container in %s", GetName());
3371  return false;
3372  }
3373 
3374  if (!h->MasterFinalize(result, fIPz, tailCut,tailMax)) {
3375  AliWarningF("Failed to finalize %s/%s", GetName(), h->GetName());
3376  return false;
3377  }
3378 
3379  Container* bgCont = GetC(result, h->GetName());
3380  if (!bgCont) {
3381  AliWarningF("%s/%s didn't put a container on output",
3382  GetName(), h->GetName());
3383  return false;
3384  }
3385  const char* sub = h->GetName();
3386  TString shrt(Form("%c", sub[0]));
3387  shrt.ToUpper();
3388  if (genCont) shrt.Append("'");
3389 
3390  TH2* background = 0;
3391  TH2* etaIPzScale = 0;
3392  if (h->GetMask() == AliAODTracklet::kInjection) {
3393  // Get the tail ratio per eta, ipZ
3394  // We get the integral of the measured eta,IPz,Delta dist
3395  // We get the integral of the injected eta,IPz,Delta dist
3396  // Then we take the ratio
3397  etaIPzScale = CopyH2(measCont, "etaIPzDeltaTail",
3398  "etaIPzScale");
3399  etaIPzScale->Divide(GetH2(bgCont, "etaIPzDeltaTail"));
3400  etaIPzScale->SetZTitle("k_{#eta,IP_{#it{z}}}");
3401  etaIPzScale->SetTitle(Form("k_{%s,#eta,IP_{#it{z}}}",shrt.Data()));
3402  bgCont->Add(etaIPzScale);
3403 
3404  // Get the tail ratio per eta
3405  // We get the integral of the measured eta,Delta dist
3406  // We get the integral of the injected eta,Detla dist
3407  // Then we take the ratio
3408  TH1* etaScale = CopyH1(measCont, "etaDeltaTail",
3409  "etaScale");
3410  etaScale->Divide(GetH1(bgCont, "etaDeltaTail"));
3411  etaScale->SetYTitle("k_{#eta}");
3412  etaScale->SetTitle(Form("k_{%s,#eta}", shrt.Data()));
3413  bgCont->Add(etaScale);
3414 
3415  // Get the integrated tail ratio.
3416  // We get the integrated tail of the measured delta dist
3417  // We get the integrated tail of the ianjection delta dist
3418  // We then get the ratio.
3419  Double_t measIntg = GetD(measCont, "deltaTail", -1);
3420  Double_t measIntE = GetD(measCont, "deltaTailError", -1);
3421  Double_t bgIntg = GetD(bgCont, "deltaTail", -1);
3422  Double_t bgIntE = GetD(bgCont, "deltaTailError", -1);
3423  Double_t scaleE = 0;
3424  Double_t scale = RatioE(measIntg, measIntE, bgIntg, bgIntE, scaleE);
3425  bgCont->Add(new TParameter<double>("scale", scale));
3426  bgCont->Add(new TParameter<double>("scaleError", scaleE));
3427 
3428  // Get the fully differential Delta distribution and scale by the
3429  // eta,IPz scalar.
3430  TH3* scaledEtaDeltaIPz = ScaleDelta(CopyH3(bgCont, "etaDeltaIPz",
3431  "scaleEtaDeltaIPz"),
3432  etaIPzScale);
3433  scaledEtaDeltaIPz->SetTitle(Form("%5.3f#times%s",
3434  scale, scaledEtaDeltaIPz->GetTitle()));
3435  scaledEtaDeltaIPz->SetDirectory(0);
3436  scaledEtaDeltaIPz->SetYTitle("k#timesd^{3}#it{N}/"
3437  "(d#Deltad#etadIP_{#it{z}})");
3438  bgCont->Add(scaledEtaDeltaIPz);
3439 #if 0
3440  // scale by derived scalars, rather than by taking the scaled full
3441  // distribution.
3442  TH2* scaledEtaDelta = CopyH2(bgCont, "etaDelta", "scaledEtaDelta");
3443  scaledEtaDelta->SetTitle(scaledEtaDeltaIPz->GetTitle());
3444  scaledEtaDelta->SetZTitle("k#timesd^{2}#it{N}/(d#Deltad#eta)");
3445  Scale(scaledEtaDelta, etaScale);
3446  bgCont->Add(scaledEtaDelta);
3447 
3448  TH1* scaledDelta = CopyH1(bgCont, "delta", "scaledDelta");
3449  scaledDelta->SetTitle(scaledEtaDeltaIPz->GetTitle());
3450  scaledDelta->SetYTitle("k#timesd#it{N}/d#Delta");
3451  Scale(scaledDelta,scale,scaleE);
3452  bgCont->Add(scaledDelta);
3453 #else
3454  // Make 2D projection to eta,Delta
3455  TH2* scaledEtaDelta = ProjectEtaDelta(scaledEtaDeltaIPz);
3456  scaledEtaDelta->SetName("scaledEtaDelta");
3457  scaledEtaDelta->SetTitle(scaledEtaDeltaIPz->GetTitle());
3458  scaledEtaDelta->SetYTitle("k#timesd^{2}#it{N}/(d#Deltad#eta)");
3459  bgCont->Add(scaledEtaDelta);
3460 
3461  // Make projection of delta
3462  TH1* scaledDelta = ProjectDelta(scaledEtaDelta);
3463  scaledDelta->SetName("scaledDelta");
3464  scaledDelta->SetTitle(scaledEtaDeltaIPz->GetTitle());
3465  scaledDelta->SetYTitle("k#timesd#it{N}/d#Delta");
3466  bgCont->Add(scaledDelta);
3467 #endif
3468  // Make background scaled by full tail
3469  background = CopyH2(bgCont, "etaIPz", "background");
3470  if (!background) AliWarningF("Didn't get background in %s", sub);
3471  else background->Multiply(etaIPzScale);
3472  // else Scale(background, scale, scaleE);
3473  }
3474  else {
3475  // For non-injection sets (i.e., combinatorics) we cannot form the
3476  // scaled background until we have the real data, so we calculate
3477  // beta instead.
3478  background = CopyH2(bgCont, "etaIPz", "background");
3479  TH2* beta = CopyH2(bgCont, "etaIPz", "beta");
3480  if (!background || !beta)
3481  AliWarningF("Didn't get background or beta in %s", sub);
3482  else {
3483  beta->Divide(GetH1(measCont, "etaIPz"));
3484  beta->SetTitle(Form("#beta_{%s}", shrt.Data()));
3485  bgCont->Add(beta);
3486  }
3487  }
3488  if (!background) {
3489  AliWarningF("Didn't get background in %s", sub);
3490  return false;
3491  }
3492  background->SetTitle(Form("#it{B}_{%s}", shrt.Data()));
3493  bgCont->Add(background);
3494 
3495  TH2* signal = CopyH2(measCont, "etaIPz", "signal");
3496  if (!signal) {
3497  AliWarningF("Didn't get signal in %s", sub);
3498  return false;
3499  }
3500  else {
3501  signal->SetTitle(Form("#it{S}_{%s}", shrt.Data()));
3502  signal->Add(background,-1);
3503  // Zero small bins
3504  for (Int_t i = 1; i <= signal->GetNbinsX(); i++) {
3505  for (Int_t j = 1; j <= signal->GetNbinsX(); j++) {
3506  if (signal->GetBinContent(i,j)<1e-6) {
3507  signal->SetBinContent(i,j,0);
3508  signal->SetBinError (i,j,0);
3509  }
3510  }
3511  }
3512  CopyAttr(background, signal);
3513  bgCont->Add(signal);
3514  }
3515 
3516  TH1* alpha = 0;
3517  if (genCont) {
3518  alpha = CopyH2(genCont, "etaIPz", "alpha");
3519  if (alpha && signal) {
3520  alpha->Divide(signal);
3521  alpha->SetTitle(Form("#alpha_{%s}", shrt.Data()));
3522  CopyAttr(signal, alpha);
3523  bgCont->Add(alpha);
3524  }
3525  }
3526 
3527  return true;
3528 }
3529 
3530 //====================================================================
3533  const char* weights,
3534  const char* sumFile,
3535  const char* resFile)
3536 {
3537  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
3538  if (!mgr) {
3539  ::Error("Create","No analysis manager to connect to.");
3540  return 0;
3541  }
3542  AliTrackletAODdNdeta* ret = 0;
3543  if (mc) {
3544  if (weights && weights[0] != '\0') {
3546  new AliTrackletAODWeightedMCdNdeta("MidRapidityMC");
3547  TUrl wurl(weights);
3548  TFile* wfile = TFile::Open(wurl.GetFile());
3549  if (!wfile) {
3550  ::Warning("Create", "Failed to open weights file: %s",
3551  wurl.GetUrl());
3552  return 0;
3553  }
3554  TString wnam(wurl.GetAnchor());
3555  if (wnam.IsNull()) wnam = "weights";
3556  TObject* wobj = wfile->Get(wnam);
3557  if (!wobj) {
3558  ::Warning("Create", "Failed to get weights %s from file %s",
3559  wnam.Data(), wfile->GetName());
3560  return 0;
3561  }
3562  if (!wobj->IsA()->InheritsFrom(AliTrackletBaseWeights::Class())) {
3563  ::Warning("Create", "Object %s from file %s not an "
3564  "AliTrackletBaseWeights but a %s",
3565  wnam.Data(), wfile->GetName(), wobj->ClassName());
3566  return 0;
3567  }
3568  wret->SetWeights(static_cast<AliTrackletBaseWeights*>(wobj));
3569  ret = wret;
3570  }
3571  else
3572  ret = new AliTrackletAODMCdNdeta("MidRapidityMC");
3573  }
3574  else ret = new AliTrackletAODdNdeta("MidRapidity");
3575  if (ret) ret->Connect();
3576 
3577  return ret;
3578 
3579 }
3580 
3581 
3582 
3583 
3584 //____________________________________________________________________
virtual Bool_t ProcessTracklet(AliAODTracklet *tracklet, Double_t ipz, UShort_t signal, Double_t weight)=0
Int_t color[]
static TProfile * GetP1(Container *parent, const char *name, Bool_t verb=true)
virtual void SetWeights(AliTrackletBaseWeights *w)
Int_t pdg
Histos & operator=(const Histos &)
static const TAxis & PdgAxis()
static void PdgAttr(Int_t pdg, TString &nme, Color_t &c, Style_t &s)
static TH2 * Make2D(Container *c, const TString &name, const TString &title, Color_t color, Style_t style, const TAxis &xAxis, const TAxis &yAxis)
virtual void SetWeightCalc(UChar_t mode=0)
void Print(Option_t *option="") const
Real_t GetEta() const
double Double_t
Definition: External.C:58
Bool_t Connect(const char *sumFile=0, const char *resFile=0)
const AliVVertex * CheckIP(const AliVVertex *ip, Double_t maxDispersion=0.04, Double_t maxZError=0.25)
virtual Bool_t FinalizeInit(Container *parent)=0
virtual void SetWeightMask(UChar_t mask=0xFF)
AliTrackletAODWeightedMCdNdeta & operator=(const AliTrackletAODWeightedMCdNdeta &o)
Bool_t MasterFinalize(Container *parent, TH1 *ipz, Double_t tailCut, Double_t tailMax)
Bool_t FinalizeInit(Container *parent)
Definition: External.C:244
static void PrintAxis(const TAxis &axis, Int_t nSig=2, const char *alt=0)
Double_t FindMultCentrality(AliVEvent *event, Int_t &nTracklets)
static void SetAxis(TAxis &axis, Int_t n, Double_t *borders)
virtual Double_t LookupWeight(AliAODTracklet *tracklet, Double_t cent, Double_t ipz)
Bool_t ProjectEtaDeltaPdg(Container *result, Int_t nEvents, Double_t tailCut, Double_t tailMax)
void SetCalc(UChar_t mode=kProduct)
Double_t FindCompatCentrality(AliVEvent *event)
virtual void SetWeightInverse(Bool_t inv)
static TProfile2D * Make2P(Container *c, const TString &name, const TString &title, Color_t color, Style_t style, const TAxis &xAxis, const TAxis &yAxis)
static TH1 * GetH1(Container *parent, const char *name, Bool_t verb=true)
static const char * CentName(Double_t c1, Double_t c2)
void SetEtaAxis(Int_t n, Double_t max)
Bool_t ProjectEtaDeltaPdgPart(Container *result, Int_t nEvents, Double_t tailDelta, Double_t tailMax, const TString &pre, const TString &tit)
static TH3 * ScaleDelta(TH3 *h, TH2 *scale)
Bool_t WorkerInit(Container *parent, const TAxis &etaAxis, const TAxis &ipzAxis, const TAxis &deltaAxis)
void SetInverse(Bool_t inv)
TCanvas * c
Definition: TestFitELoss.C:172
virtual void SetWeightVeto(UChar_t veto=0xFF)
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 SetMaxDelta(Double_t x=25)
ClassDef(AliTrackletAODdNdeta, 1)
static TH2 * ScaleToIPz(TH2 *h, TH1 *ipZ, Bool_t full=false)
Int_t nCentBins
static TH2 * ProjectEtaDelta(TH3 *h)
void SetCentralityAxis(const TString &spec)
static TH1 * AverageOverIPz(TH2 *h, const char *name, UShort_t mode, TH1 *ipz, TH2 *mask=0, Bool_t verb=true)
static TH1 * ProjectDelta(TH2 *h)
void SetWeightVeto(UChar_t veto=0x0)
static Double_t GetD(Container *parent, const char *name, Double_t def=-1, Bool_t verb=true)
virtual Real_t GetParentPt(Bool_t second=false) const
void SetTailDelta(Double_t x=5)
static TH3 * Make3D(Container *c, const TString &name, const TString &title, Color_t color, Style_t style, const TAxis &xAxis, const TAxis &yAxis, const TAxis &zAxis)
static TH1 * Scale(TH1 *h, Double_t x, Double_t xe)
Definition: External.C:92
void SetMask(UChar_t mask)
Bool_t IsCombinatorics() const
CentBin & operator=(const CentBin &)
virtual Bool_t WorkerInit(Container *parent, const TAxis &etaAxis, const TAxis &ipzAxis, const TAxis &deltaAxis)
Bool_t CheckEvent(Double_t &cent, const AliVVertex *&ip, TClonesArray *&tracklets)
void SetIPzAxis(Int_t n, Double_t max)
Utilities for midrapidity analysis.
void SetAbsMinCent(Double_t x=-1)
Tracklet AOD object.
CentBin & operator=(const CentBin &)
void Print(Option_t *option="") const
TClonesArray * FindTracklets(AliVEvent *event)
Bool_t MasterFinalize(Container *parent, TH1 *ipz, Double_t tailCut, Double_t tailMax)
ClassDef(AliTrackletAODMCdNdeta, 1)
void SetTailMaximum(Double_t x=-1)
void SetEtaAxis(Int_t n, Double_t min, Double_t max)
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
Bool_t InitCentBins(Container *existing)
float Float_t
Definition: External.C:68
virtual void SetDebug(UShort_t lvl)
Real_t GetDPhi() const
static Double_t Integrate(TH1 *h, Double_t min, Double_t max, Double_t &err)
virtual Bool_t MasterFinalize(Container *parent, TH1 *ipz, Double_t tailCut, Double_t tailMax)=0
void SetCentralityMethod(const TString &name)
Bool_t FinalizeInit(Container *parent)
Real_t GetPhi() const
Encode simulation weights for 2nd pass.
const AliVVertex * FindIP(AliVEvent *event, Double_t maxDispersion=0.04, Double_t maxZError=0.25)
Definition: External.C:212
const AliVVertex * FindSimpleIP(AliVEvent *event, Double_t maxDispersion=0.04, Double_t maxZError=0.25)
void SetDPhiShift(Double_t x=0.0045)
AliTrackletAODMCdNdeta(const AliTrackletAODdNdeta &o)
Double_t FindCentrality(AliVEvent *event, Int_t &nTracklets)
virtual AliTrackletAODdNdeta::CentBin * MakeCentBin(Float_t c1, Float_t c2)
void SetShiftedDPhiCut(Double_t x=-1)
void Print(Option_t *option="") const
void SetIPzAxis(Int_t n, Double_t min, Double_t max)
const AliVVertex * FindRealIP(AliVEvent *event, Double_t maxDispersion=0.04, Double_t maxZError=0.25)
Bool_t IsMeasured() const
static TH2 * CopyH2(Container *parent, const char *name, const char *newName=0, Bool_t verb=true)
virtual Short_t GetParentPdg(Bool_t second=false) const
AliTrackletAODWeightedMCdNdeta(const char *name)
virtual void Print(Option_t *option="") const
Bool_t EstimateBackground(Container *result, Container *measCont, Container *genCont, Histos *h, Double_t tailCut, Double_t tailMax)
void SetDeltaCut(Double_t x=1.5)
Int_t mode
Definition: anaM.C:40
static TH1 * CopyH1(Container *parent, const char *name, const char *newName=0, Bool_t verb=true)
virtual const char * GetBranchName() const
void SetVeto(UChar_t veto)
static TH3 * GetH3(Container *parent, const char *name, Bool_t verb=true)
const char * GetName() const
static Double_t RatioE(Double_t n, Double_t en, Double_t d, Double_t ed, Double_t &er)
static TH2 * GetH2(Container *parent, const char *name, Bool_t verb=true)
static AliTrackletAODdNdeta * Create(Bool_t mc=false, const char *weights=0, const char *sumFile=0, const char *resFile=0)
UChar_t GetFlags() const
Histos * MakeHistos(const char *name, Color_t color, Style_t style, UShort_t mask, UShort_t veto)
Bool_t IsNeutral() const
AliTrackletAODMCdNdeta(const char *name)
Bool_t MasterFinalize(Container *results)
void SetPhiAxis(Int_t n, Double_t min, Double_t max)
Bool_t IsGenerated() const
Bool_t ProjectEtaPdgIPz(Container *result, TH1 *ipz, const TString &shn)
ClassDef(AliTrackletAODWeightedMCdNdeta, 1)
Float_t nEvents[nProd]
void SetAttr(Color_t c, Style_t m)
Definition: External.C:220
void SetEtaAxis(const TString &spec)
virtual Double_t LookupWeight(AliAODTracklet *tracklet, Double_t cent, Double_t ipz)
virtual void Print(Option_t *option="") const
static Container * GetC(Container *parent, const char *name, Bool_t verb=true)
static Int_t PdgBin(Int_t pdg)
AliTrackletAODWeightedMCdNdeta(const AliTrackletAODWeightedMCdNdeta &o)
Bool_t ProjectEtaPdg(Container *result, Int_t nEvents)
Bool_t ProcessTracklet(AliAODTracklet *tracklet, Double_t ipz, UShort_t signal, Double_t weight)
Bool_t Accept(Double_t cent, Double_t ipz)
unsigned short UShort_t
Definition: External.C:28
void Print(Option_t *option="") const
virtual UShort_t CheckTracklet(AliAODTracklet *tracklet) const
void SetIPzAxis(const TString &spec)
const char Option_t
Definition: External.C:48
virtual Bool_t MasterFinalize(Container *results)
Real_t GetDelta() const
void SetWeightMask(UChar_t mask=0xFF)
AliTrackletBaseWeights * fWeights
static TH3 * CopyH3(Container *parent, const char *name, const char *newName=0, Bool_t verb=true)
Histos(const char *name="", Color_t color=kBlack, Style_t style=1, UChar_t mask=0, UChar_t veto=0)
bool Bool_t
Definition: External.C:53
static TObject * CloneAndAdd(Container *c, TObject *o)
virtual void SetDebug(UShort_t lvl)
virtual const char * GetBranchName() const
Bool_t WorkerInit(Container *parent, const TAxis &etaAxis, const TAxis &ipzAxis, const TAxis &deltaAxis)
static TH1 * Make1D(Container *c, const TString &name, const TString &title, Color_t color, Style_t style, const TAxis &xAxis)
virtual CentBin * MakeCentBin(Float_t c1, Float_t c2)
static TProfile * Make1P(Container *c, const TString &name, const TString &title, Color_t color, Style_t style, const TAxis &xAxis)
void SetPhiAxis(Int_t n, Double_t max)
static void FixAxis(TAxis &axis, const char *title=0)
A simplified AOD header.
virtual CentBin * InitCentBin(Int_t bin, Float_t c1, Float_t c2, Container *existing, const TAxis &deltaAxis)
Definition: External.C:196
void SetCentralityAxis(Int_t n, Double_t *bins)
void SetWeights(AliTrackletBaseWeights *w)
AliTrackletAODdNdeta & operator=(const AliTrackletAODdNdeta &o)
Bool_t ProcessTracklet(AliAODTracklet *tracklet, Double_t ipz, UShort_t signal, Double_t weight)
void ProcessEvent(Double_t cent, const AliVVertex *ip, TClonesArray *tracklets)
static TProfile * GetP(Container *parent, const char *name, Bool_t verb=true)
static void CopyAttr(const TH1 *src, TH1 *dest)