AliPhysics  fc9925d (fc9925d)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliTrackletdNdetaTask.C
Go to the documentation of this file.
1 #ifndef ALITRACKLETDNDETATASK_H
2 #define ALITRACKLETDNDETATASK_H
3 #include <AliAnalysisTaskSE.h>
4 #include <AliTrackletdNdetaUtils.C>
5 #ifndef __CINT__
6 #include <TGeoGlobalMagField.h>
7 #include <TProfile.h>
8 #include <AliVEvent.h>
9 #include <AliVVertex.h>
10 #include <AliVertex.h>
11 #include <AliVMultiplicity.h>
12 #include <AliMultiplicity.h>
13 #include <AliAnalysisManager.h>
14 #include <AliVEventHandler.h>
15 #include <AliESDInputHandlerRP.h>
16 #include <AliMultSelection.h>
17 #include <AliITSMultRecBg.h>
18 #include <AliCDBManager.h>
19 #include <AliCDBEntry.h>
20 #include <AliCDBPath.h>
21 #include <AliCDBId.h>
22 #include <AliGeomManager.h>
23 #include <TGeoManager.h>
24 #else
25 class AliVEvent;
26 class AliVVertex;
27 class AliMultiplicity;
28 class AliAnalysisDataContainer;
29 class AliITSMultRecBg;
30 class AliMultSelection; // Auto-load
31 class TGeoGlobalMagField; // Auto-load
32 class TGeoManager; // Auto-load
33 class AliGeomManager; // Auto-load
34 class AliCDBManager; // Auto-load
35 #endif
36 
37 namespace {
38  Int_t lDebug = 0;
39 }
45 class AliTrackletdNdetaTask : public AliAnalysisTaskSE,
46  public AliTrackletdNdetaUtils
47 {
48 public:
49  //==================================================================
53  enum {
54  kNormal = 0x1,
55  kInjection = 0x2,
56  kRotation = 0x4
57  };
58  enum {
59  kAll=1, // Count all events
60  kEvent, // Have event
61  kMultiplicity, // Have multiplicity
62  kClusters, // Have clusters
63  kTrigger, // Have trigger
64  kIP, // Have IP
65  kCentrality, // Have Centrality
66  kReconstructed, // Have reconstructed
67  kInjected, // Have injected
68  kRotated, // Have rotated
69  kCompleted // Have completed
70  };
71  //==================================================================
75  struct SubBase : public TObject
76  {
77  typedef AliTrackletdNdetaTask::Container Container;
79  Container* fContainer;
83  SubBase() : fContainer(0) {}
87  SubBase(const SubBase& o) : TObject(o), fContainer(0) {}
91  virtual ~SubBase() {}
97  SubBase& operator=(SubBase& o) { return *this; }
103  virtual const char* Name() const = 0;
109  virtual Container* CreateContainer(Container& parent);
120  virtual void WorkerInit(Container& parent,
121  const TAxis& etaAxis,
122  const TAxis& ipzAxis,
123  const TAxis& deltaAxis,
124  const TAxis& dThetaAxis,
125  const TAxis& dPhiAxis) = 0;
131  virtual void FinalizeInit(Container* parent) = 0;
142  virtual Container* MasterFinalize(Container* parent,
143  TH1* ipz,
144  Double_t deltaCut,
145  Double_t deltaTail) = 0;
146 
147  ClassDef(SubBase,1); // Base class for sub structures
148  };
149  //==================================================================
150  struct HistoSet : public SubBase
151  {
152  TString fName;
153  TH2* fEtaVsDelta; // Eta vs delta
154  TH2* fDPhiVsdThetaX; // Normalized dphi vs dtheta*x
155  TH2* fDPhiVsdTheta; // dphi vs dtheta
156  TH1* fDelta; // Distance
157  TH2* fEtaVsIPz; // Eta vs IPz
158  Color_t fColor; // Color for this histogram set
159  Style_t fStyle;
163  HistoSet(const char* name="",
164  Color_t color=kBlack,
165  Style_t style=20)
166  : SubBase(),
167  fName(name),
168  fEtaVsDelta(0),
169  fDPhiVsdTheta(0),
170  fDPhiVsdThetaX(0),
171  fDelta(0),
172  fEtaVsIPz(0),
173  fColor(color),
174  fStyle(style)
175  {}
176  HistoSet(const HistoSet& o)
177  : SubBase(o),
178  fName(o.fName),
179  fEtaVsDelta(0),
180  fDPhiVsdTheta(0),
181  fDPhiVsdThetaX(0),
182  fDelta(0),
183  fEtaVsIPz(0),
184  fColor(o.fColor),
185  fStyle(o.fStyle)
186  {}
192  HistoSet& operator=(const HistoSet&) { return *this; }
196  virtual ~HistoSet() {}
202  virtual const char* Name() const { return fName.Data(); }
213  virtual void WorkerInit(Container& parent,
214  const TAxis& etaAxis,
215  const TAxis& ipzAxis,
216  const TAxis& deltaAxis,
217  const TAxis& dThetaAxis,
218  const TAxis& dPhiAxis);
224  virtual void FinalizeInit(Container* parent);
235  virtual Container* MasterFinalize(Container* parent,
236  TH1* ipz,
237  Double_t deltaCut,
238  Double_t deltaTail);
251  virtual void Fill(Double_t ipZ,
252  Double_t eta,
253  Double_t dphi,
254  Double_t dphis,
255  Double_t dtheta,
256  Double_t dthetax,
257  Double_t delta,
258  Double_t weight=1);
259  ClassDef(HistoSet,1); // Class to hold histograms
260  };
261  //==================================================================
265  struct CentBin : public SubBase
266  {
276  CentBin(Float_t cmin,
277  Float_t cmax,
278  UShort_t recFlags,
279  Int_t colOff=1,
280  Int_t styOff=0);
284  virtual ~CentBin() {}
291  virtual const char* Name() const;
303  virtual void WorkerInit(Container& parent,
304  const TAxis& etaAxis,
305  const TAxis& ipzAxis,
306  const TAxis& deltaAxis,
307  const TAxis& dThetaAxis,
308  const TAxis& dPhiAxis);
314  virtual void FinalizeInit(Container* parent);
324  virtual Container* EstimateBackground(Container* dataCon,
325  HistoSet* set,
326  Container* result,
327  Double_t deltaCut,
328  Double_t deltaTail);
339  virtual Container* MasterFinalize(Container* parent,
340  TH1* ipz,
341  Double_t deltaCut,
342  Double_t deltaTail);
354  virtual Bool_t Accept(Float_t value) const
355  {
356  return (fMin <= value && value < fMax);
357  }
363  virtual void FillIPz(Double_t ipZ) { fIPz->Fill(ipZ); }
381  virtual Bool_t FillSet(UShort_t set,
382  Bool_t isSignal,
383  Double_t ipZ,
384  Double_t eta,
385  Double_t dPhi,
386  Double_t dPhiS,
387  Double_t dTheta,
388  Double_t dThetaX,
389  Double_t delta,
390  Double_t weight);
391  /* @} */
392  protected:
401  CentBin()
402  : SubBase(),
403  fMin(0),
404  fMax(0),
405  fHistoSets(0),
406  fDataSet(0),
407  fInjSet(0),
408  fRotSet(0),
409  fIPz(0)
410  {}
414  CentBin(const CentBin& o)
415  : SubBase(o),
416  fMin(o.fMin),
417  fMax(o.fMax),
418  fHistoSets(0),
419  fDataSet(0),
420  fInjSet(0),
421  fRotSet(0),
422  fIPz(0)
423  {
424  }
430  CentBin& operator=(const CentBin&) { return *this; }
431  /* @} */
432 
437  Float_t fMin; // Least value
438  Float_t fMax; // Largest value
439  TList* fHistoSets; // List of histogram sets
440  HistoSet* fDataSet; // Normal reconstruction histograms
441  HistoSet* fInjSet; // Injected histograms
442  HistoSet* fRotSet; // Rotation histograms
443  TH1* fIPz; // Vertex distribution
444  /* @} */
445  ClassDef(CentBin,1); // Class to centrality dependent stuff
446  };
447  //==================================================================
457  AliTrackletdNdetaTask(const char* name);
461  AliTrackletdNdetaTask()
462  : AliAnalysisTaskSE(),
463  AliTrackletdNdetaUtils(),
464  fContainer(0),
465  fCentBins(0),
466  fCentMethod("V0M"),
467  fCentAxis(1,0,100),
468  fIPzAxis(20,-10,10),
469  fEtaAxis(80,-2,2),
470  fPhiAxis(200,0,TMath::TwoPi()),
471  fCent(0),
472  fIPz(0),
473  fEtaPhi(0),
474  fUsedClusters0(0),
475  fUsedClusters1(0),
476  fAllClusters0(0),
477  fAllClusters1(0),
478  fStatus(0),
479  fRecMode(kNormal|kInjection),
480  fScaleDTheta(false),
481  fDPhiShift(0.0045),
482  fShiftedDPhiCut(-1),
483  fScaledDThetaCut(-1),
484  fDeltaCut(1.5),
485  fMaxDelta(25),
486  fTailDelta(5),
487  fDThetaWindow(0.025),
488  fDPhiWindow(0.06),
489  fPhiOverlapCut(0.005),
490  fZEtaOverlapCut(0.05),
491  fPhiRotation(TMath::Pi())
492  {}
496  virtual ~AliTrackletdNdetaTask() {}
497  //__________________________________________________________________
506  void UserCreateOutputObjects();
510  void UserExec(Option_t*);
515  void FinishTaskOutput() { WorkerFinalize(); }
520  void Terminate(Option_t*);
529  virtual Bool_t Connect(const char* sumFile=0,
530  const char* resFile=0);
535  void SetReconstructionMode(const TString& mode)
536  {
537  UShort_t flags = 0;
538  TString m(mode); m.ToLower();
539  if (m.Contains("nor")) flags |= kNormal;
540  if (m.Contains("inj")) flags |= kInjection;
541  if (m.Contains("rot")) flags |= kRotation;
542  SetReconstructionMode(flags);
543  }
549  void SetReconstructionMode(UShort_t flags) { fRecMode = flags; }
555  void SetCentralityMethod(const TString& name) { fCentMethod = name; }
562  void SetCentralityAxis(Int_t n, Double_t* bins)
563  {
564  SetAxis(fCentAxis,n,bins);
565  }
571  void SetCentralityAxis(const TString& spec)
572  {
573  SetAxis(fCentAxis, spec, "-:,");
574  }
582  void SetIPzAxis(Int_t n, Double_t min, Double_t max)
583  {
584  SetAxis(fIPzAxis, n, min, max);
585  }
592  void SetIPzAxis(Int_t n, Double_t max)
593  {
594  SetAxis(fIPzAxis, n, max);
595  }
596  void SetIPzAxis(const TString& spec)
597  {
598  SetAxis(fIPzAxis, spec);
599  }
607  void SetPhiAxis(Int_t n, Double_t min, Double_t max)
608  {
609  SetAxis(fPhiAxis, n, min, max);
610  }
617  void SetPhiAxis(Int_t n, Double_t max)
618  {
619  SetAxis(fPhiAxis, n, max);
620  }
628  void SetEtaAxis(Int_t n, Double_t min, Double_t max)
629  {
630  SetAxis(fEtaAxis, n, min, max);
631  }
638  void SetEtaAxis(Int_t n, Double_t max)
639  {
640  SetAxis(fEtaAxis, n, max);
641  }
642  void SetEtaAxis(const TString& spec)
643  {
644  SetAxis(fEtaAxis, spec);
645  }
651  void SetScaleDTheta(Bool_t x=false) { fScaleDTheta = x; }
657  void SetDPhiShift(Double_t x=0.0045) { fDPhiShift = x; }
663  void SetShiftedDPhiCut(Double_t x=-1) { fShiftedDPhiCut = x; }
669  void SetScaleDThetaCut(Double_t x=-1) { fScaledDThetaCut = x; }
675  void SetDeltaCut(Double_t x=1.5) { fDeltaCut = x; }
681  void SetMaxDelta(Double_t x=25) { fMaxDelta = x; }
687  void SetTailDelta(Double_t x=5) { fTailDelta = x; }
693  void SetDThetaWindow(Double_t x=0.025) { fDThetaWindow = x; }
699  void SetDPhiWindow(Double_t x=0.06) { fDPhiWindow = x; }
705  void SetPhiOverlapCut(Double_t x=0.005) { fPhiOverlapCut = x; }
711  void SetZEtaOverlapCut(Double_t x=0.05) { fZEtaOverlapCut = x; }
717  void SetPhiRotation(Double_t x=TMath::Pi()) { fPhiRotation = x; }
718  /* @} */
723  virtual void Print(Option_t*) const;
724 protected:
725  //__________________________________________________________________
735  AliTrackletdNdetaTask(const AliTrackletdNdetaTask& o)
736  : AliAnalysisTaskSE(o),
737  AliTrackletdNdetaUtils(o),
738  fContainer(0),
739  fCentBins(0),
740  fCentMethod("V0M"),
741  fCentAxis(1,0,100),
742  fIPzAxis(20,-10,10),
743  fEtaAxis(80,-2,2),
744  fPhiAxis(200,0,TMath::TwoPi()),
745  fCent(0),
746  fIPz(0),
747  fEtaPhi(0),
748  fUsedClusters0(0),
749  fUsedClusters1(0),
750  fAllClusters0(0),
751  fAllClusters1(0),
752  fStatus(0),
753  fRecMode(kNormal|kInjection),
754  fScaleDTheta(false),
755  fDPhiShift(0.0045),
756  fShiftedDPhiCut(-1),
757  fScaledDThetaCut(-1),
758  fDeltaCut(1.5),
759  fMaxDelta(25),
760  fTailDelta(5),
761  fDThetaWindow(0.025),
762  fDPhiWindow(0),
763  fPhiOverlapCut(0.005),
764  fZEtaOverlapCut(0.05),
765  fPhiRotation(TMath::Pi())
766  {}
774  AliTrackletdNdetaTask& operator=(const AliTrackletdNdetaTask& o)
775  {
776  if (&o == this) return *this;
777  return *this;
778  }
780  //__________________________________________________________________
787  virtual void InitCentBins(Container* existing);
793  virtual Bool_t InitCDB();
799  virtual Bool_t InitGeometry();
805  virtual Int_t GetCDBReferenceRun() const { return 137161; }
811  virtual const char* GetCDBReferenceURL() const
812  {
813  return "alien://Folder=/alice/data/2010/OCDB";
814  }
819  virtual void WorkerInit();
824  virtual void WorkerFinalize();
830  virtual void MasterFinalize(Container* results);
831  /* @} */
841  virtual TTree* FindClusters(Bool_t needed);
850  virtual Double_t FindCentrality(AliVEvent* event);
860  virtual const AliVVertex* FindIP(AliVEvent* event,
861  Double_t maxDispersion=0.04,
862  Double_t maxZError=0.25);
868  virtual Bool_t FindTrigger();
877  virtual Double_t TrackletWeight(Int_t trackletNumber,
878  const AliMultiplicity* mult) const
879  {
880  return 1;
881  }
893  virtual Bool_t ProcessTracklets(TList& toRun,
894  AliITSMultRecBg* reco,
895  Double_t cent,
896  Double_t ipZ,
897  AliMultiplicity* mult);
905  UShort_t FindMode(AliITSMultRecBg* reco) const;
923  virtual void FillBins(TList& toRun,
924  AliITSMultRecBg* reco,
925  AliMultiplicity* mult,
926  Int_t trackletNumber,
927  Bool_t isSignal,
928  Double_t ipz,
929  Double_t eta,
930  Double_t dPhi,
931  Double_t dPhiS,
932  Double_t dTheta,
933  Double_t dThetaX,
934  Double_t delta,
935  Double_t weight) const;
945  virtual Bool_t Reconstruct(TList& toRun,
946  UShort_t mode,
947  TTree* clusters,
948  Double_t cent,
949  const AliVVertex* ip);
956  virtual void FillClusters(AliITSMultRecBg* reco);
964  virtual void FillClusters(AliMultiplicity* mult, Double_t ipZ);
976  virtual Bool_t CheckEvent(Double_t& cent,
977  const AliVVertex*& ip,
978  AliMultiplicity*& mult,
979  TTree*& clusters);
989  virtual void ProcessEvent(Double_t cent,
990  const AliVVertex* ip,
991  AliMultiplicity* mult,
992  TTree* clusters);
1002  virtual CentBin* MakeCentBin(Float_t c1, Float_t c2, UShort_t recFlags)
1003  {
1004  return new CentBin(c1, c2, recFlags);
1005  }
1006  /* @} */
1007  //__________________________________________________________________
1012  Container* fContainer;
1013  Container* fCentBins;
1014  TString fCentMethod;
1015  TAxis fCentAxis;
1016  TAxis fIPzAxis;
1017  TAxis fEtaAxis;
1018  TAxis fPhiAxis;
1019  TH1* fCent;
1020  TH1* fIPz;
1021  TH2* fEtaPhi;
1022  TH2* fUsedClusters0;
1023  TH2* fUsedClusters1;
1024  TH2* fAllClusters0;
1025  TH2* fAllClusters1;
1026  TH1* fStatus;
1027  UShort_t fRecMode;
1028  Bool_t fScaleDTheta;
1029  Double_t fMaxDelta;
1030  Double_t fTailDelta;
1031  Double_t fDPhiShift;
1032  Double_t fShiftedDPhiCut;
1033  Double_t fScaledDThetaCut;
1034  Double_t fDeltaCut;
1035  Double_t fDThetaWindow;
1036  Double_t fDPhiWindow;
1037  Double_t fPhiOverlapCut;
1038  Double_t fZEtaOverlapCut;
1039  Double_t fPhiRotation;
1040  /* @} */
1041 
1042  ClassDef(AliTrackletdNdetaTask,1); // Base class for dN/deta u/tracklets
1043 };
1044 
1045 
1046 //====================================================================
1047 // Sub base member functions
1048 //____________________________________________________________________
1049 AliTrackletdNdetaTask::Container*
1050 AliTrackletdNdetaTask::SubBase::CreateContainer(Container& parent)
1051 {
1052  Container* c = new Container;
1053  c->SetName(Name());
1054  c->SetOwner();
1055  parent.Add(c);
1056  return c;
1057 }
1058 
1059 //====================================================================
1060 // Histogram set member functions
1061 //____________________________________________________________________
1062 void AliTrackletdNdetaTask::HistoSet::WorkerInit(Container& parent,
1063  const TAxis& etaAxis,
1064  const TAxis& ipzAxis,
1065  const TAxis& deltaAxis,
1066  const TAxis& dThetaAxis,
1067  const TAxis& dPhiAxis)
1068 {
1069  fContainer = CreateContainer(parent);
1070  Container* c = fContainer;
1071 
1072  fEtaVsIPz = Make2D(*c, "etaVsIPz",
1073  Form("#eta vs IP_{#it{z}} %s",Name()),
1074  fColor,fStyle, etaAxis, ipzAxis);
1075  fEtaVsDelta = Make2D(*c, "etaVsDelta","",fColor,fStyle,etaAxis,deltaAxis);
1076  fDPhiVsdTheta = Make2D(*c, "dphiVsdtheta","",fColor,fStyle,
1077  dPhiAxis,dThetaAxis);
1078  fDPhiVsdThetaX = Make2D(*c, "dphiVsdthetaX","",fColor,fStyle,
1079  dPhiAxis,dThetaAxis);
1080  fDPhiVsdThetaX->SetXTitle("#Delta#phi-#delta#phi");
1081  fDPhiVsdThetaX->SetYTitle("#Delta#theta/sin^{2}(#theta)");
1082  fDelta = Make1D(*c, "delta", Form("#Delta %s", Name()),
1083  fColor, fStyle,deltaAxis);
1084 }
1085 //____________________________________________________________________
1086 void AliTrackletdNdetaTask::HistoSet::FinalizeInit(Container* parent)
1087 {
1088  fContainer = GetC(parent, Name());
1089  fEtaVsDelta = GetH2(fContainer, "etaVsDelta");
1090  fDPhiVsdThetaX = GetH2(fContainer, "dphiVsdthetaX");
1091  fDPhiVsdTheta = GetH2(fContainer, "dphiVsdtheta");
1092  fEtaVsIPz = GetH2(fContainer, "etaVsIPz");
1093  fDelta = GetH1(fContainer, "delta");
1094 }
1095 //____________________________________________________________________
1096 AliTrackletdNdetaTask::Container*
1097 AliTrackletdNdetaTask::HistoSet::MasterFinalize(Container* parent,
1098  TH1* ipz,
1099  Double_t deltaCut,
1100  Double_t tailDelta)
1101 {
1102  Container* result = CreateContainer(*parent);
1103 
1104  // Get the number of events
1105  Double_t nEvents = ipz->GetEntries();
1106 
1107  // Scale each vertex range by number of events in that range
1108  TH2* etaVsIPz = ScaleToIPz(fEtaVsIPz, ipz);
1109  result->Add(etaVsIPz);
1110 
1111  // Normalize delta distribution to integral over 5 to 25
1112  TH1* delta = static_cast<TH1*>(CloneAndAdd(result, fDelta));
1113  delta->Scale(1./nEvents);
1114  Double_t maxDelta = fDelta->GetXaxis()->GetXmax();
1115  Double_t eintg;
1116  Double_t intg = Integrate(delta, tailDelta, maxDelta, eintg);
1117 
1118  result->Add(new TParameter<double>("deltaTailIntg", intg));
1119  result->Add(new TParameter<double>("deltaTailIntE", eintg));
1120 
1121  TH2* etaVsDelta = static_cast<TH2*>(CloneAndAdd(result,fEtaVsDelta));
1122  TH2* dPhiVsdThetaX = static_cast<TH2*>(CloneAndAdd(result,fDPhiVsdThetaX));
1123  TH2* dPhiVsdTheta = static_cast<TH2*>(CloneAndAdd(result,fDPhiVsdTheta));
1124  etaVsDelta ->Scale(1./nEvents);
1125  dPhiVsdThetaX->Scale(1./nEvents);
1126  dPhiVsdTheta ->Scale(1./nEvents);
1127 
1128  TH1* deltaIntg = etaVsDelta->ProjectionX("deltaTailInt");
1129  deltaIntg->SetDirectory(0);
1130  deltaIntg->Reset();
1131  deltaIntg->SetYTitle(Form("#int_{%3.1f}^{%4.1f}d#Delta",
1132  tailDelta, maxDelta));
1133  for (Int_t i = 1; i <= deltaIntg->GetNbinsX(); i++) {
1134  TH1* tmp = etaVsDelta->ProjectionY("tmp",i,i);
1135  intg = Integrate(tmp, tailDelta, maxDelta, eintg);
1136  deltaIntg->SetBinContent(i, intg);
1137  deltaIntg->SetBinError (i, eintg);
1138  delete tmp;
1139  }
1140  result->Add(deltaIntg);
1141 
1142  return result;
1143 }
1144 //____________________________________________________________________
1145 void AliTrackletdNdetaTask::HistoSet::Fill(Double_t ipZ,
1146  Double_t eta,
1147  Double_t dphi,
1148  Double_t dphis,
1149  Double_t dtheta,
1150  Double_t dthetax,
1151  Double_t delta,
1152  Double_t weight)
1153 {
1154  if (delta > fEtaVsDelta->GetYaxis()->GetXmax()) return;
1155  fEtaVsDelta ->Fill(eta, delta, weight);
1156  fDPhiVsdThetaX ->Fill(dphis, dthetax, weight);
1157  fDPhiVsdTheta ->Fill(dphi, dtheta, weight);
1158  fDelta ->Fill(delta, weight);
1159  if (ipZ > -999.) fEtaVsIPz->Fill(eta, ipZ, weight);
1160 }
1161 //====================================================================
1162 // Centrality bin member functions
1163 //____________________________________________________________________
1164 AliTrackletdNdetaTask::CentBin::CentBin(Float_t cmin,
1165  Float_t cmax,
1166  UShort_t recFlags,
1167  Int_t colOff,
1168  Int_t styOff)
1169  : SubBase(),
1170  fMin(cmin),
1171  fMax(cmax),
1172  fIPz(0),
1173  fHistoSets(0),
1174  fDataSet(0),
1175  fInjSet(0),
1176  fRotSet(0)
1177 {
1178  fHistoSets = new TList;
1179  fHistoSets->SetOwner(true);
1180 
1181  fDataSet = new HistoSet("data", kRed+colOff, 20+styOff);
1182  fHistoSets->Add(fDataSet);
1183  if (recFlags & kInjection) {
1184  fInjSet = new HistoSet("injection", kGreen+colOff, 21+styOff);
1185  fHistoSets->Add(fInjSet);
1186  }
1187  if (recFlags & kRotation) {
1188  fRotSet = new HistoSet("rotation", kBlue+colOff, 22+styOff);
1189  fHistoSets->Add(fRotSet);
1190  }
1191 }
1192 //____________________________________________________________________
1193 const char* AliTrackletdNdetaTask::CentBin::Name() const
1194 {
1195  return Form("%s%03dd%02d_%03dd%02d", "cent",
1196  Int_t(fMin), Int_t(fMin*100) % 100,
1197  Int_t(fMax), Int_t(fMax*100) % 100);
1198 }
1199 //____________________________________________________________________
1200 void AliTrackletdNdetaTask::CentBin::WorkerInit(Container& parent,
1201  const TAxis& etaAxis,
1202  const TAxis& ipzAxis,
1203  const TAxis& deltaAxis,
1204  const TAxis& dThetaAxis,
1205  const TAxis& dPhiAxis)
1206 {
1207  fContainer = CreateContainer(parent);
1208  fIPz = Make1D(*fContainer, "ipz", "", kMagenta+2, 20, ipzAxis);
1209  TIter next(fHistoSets);
1210  HistoSet* set = 0;
1211  while ((set = static_cast<HistoSet*>(next())))
1212  set->WorkerInit(*fContainer,etaAxis,ipzAxis,deltaAxis,dThetaAxis,dPhiAxis);
1213 }
1214 //____________________________________________________________________
1215 void AliTrackletdNdetaTask::CentBin::FinalizeInit(Container* parent)
1216 {
1217  fContainer = GetC(parent, Name());
1218  fIPz = GetH1(fContainer, "ipz");
1219  TIter next(fHistoSets);
1220  HistoSet* set = 0;
1221  while ((set = static_cast<HistoSet*>(next())))
1222  set->FinalizeInit(fContainer);
1223 }
1224 //____________________________________________________________________
1225 AliTrackletdNdetaTask::Container*
1226 AliTrackletdNdetaTask::CentBin::EstimateBackground(Container* dataCon,
1227  HistoSet* set,
1228  Container* result,
1229  Double_t deltaCut,
1230  Double_t deltaTail)
1231 {
1232  if (!set || !dataCon) return 0;
1233 
1234  Container* bgCon = set->MasterFinalize(result, fIPz, deltaCut, deltaTail);
1235  if (!bgCon) return 0;
1236 
1237  Bool_t isComb = set->fName.Contains("combinatorics", TString::kIgnoreCase);
1238  // Info("EstimateBackground", "%s - combinatorics? %s",
1239  // set->fName.Data(), isComb ? "yes" : "no");
1240 
1241  Double_t dataIntg = GetD(dataCon, "deltaTailIntg", -1);
1242  Double_t dataIntE = GetD(dataCon, "deltaTailIntE", -1);
1243  Double_t bgIntg = GetD(bgCon, "deltaTailIntg", -1);
1244  Double_t bgIntE = GetD(bgCon, "deltaTailIntE", -1);
1245  if (dataIntg <= 0 || bgIntg <= 0) return bgCon;
1246  if (lDebug > 2) {
1247  AliInfoF("int(measured) Delta tail: %f +/- %f", dataIntg, dataIntE);
1248  AliInfoF("int(%s) Delta tail: %f +/- %f", bgCon->GetName(),
1249  bgIntg, bgIntE);
1250  }
1251  Double_t scaleE = 0;
1252  Double_t scale = RatioE(dataIntg, dataIntE, bgIntg, bgIntE, scaleE);
1253 
1254  TH1* bgDelta = GetH1(bgCon, "delta");
1255  if (!bgDelta) return bgCon;
1256 
1257  TH1* dataDelta = GetH1(dataCon, "delta");
1258  if (!dataDelta) {
1259  AliWarningF("No delta distribution in %s", dataCon->GetName());
1260  }
1261 
1262  TH1* deltaScaled = CopyH1(bgCon, "delta", "deltaScaled");
1263  bgCon->Add(deltaScaled);
1264  deltaScaled->SetLineStyle(7);
1265  deltaScaled->SetTitle(Form("%s#times%6.3f", bgDelta->GetTitle(), scale));
1266  if (dataDelta) deltaScaled->SetMarkerStyle(dataDelta->GetMarkerStyle());
1267  else deltaScaled->SetMarkerStyle(20);
1268 
1269  if (lDebug > 2)
1270  AliInfoF("Scaling %s Delta dist by %f +/- %f",
1271  bgCon->GetName(), scale, scaleE);
1272  Scale(deltaScaled, scale, scaleE);
1273 
1274  TH2* backgroundEst = CopyH2(bgCon, "etaVsIPz", "backgroundEst");
1275  backgroundEst->SetTitle("Background");
1276  if (!isComb) Scale(backgroundEst, scale, scaleE); // ->Scale(scale);
1277  // else Info("EstimateBackground", "Combinators, no scaling of BG");
1278  bgCon->Add(backgroundEst);
1279 
1280  TH2* signalEst = CopyH2(dataCon, "etaVsIPz","signalEst");
1281  signalEst->SetMarkerStyle(backgroundEst->GetMarkerStyle());
1282  signalEst->SetMarkerColor(backgroundEst->GetMarkerColor());
1283  signalEst->SetMarkerSize (backgroundEst->GetMarkerSize());
1284  signalEst->SetLineStyle (backgroundEst->GetLineStyle());
1285  signalEst->SetLineColor (backgroundEst->GetLineColor());
1286  signalEst->SetLineWidth (backgroundEst->GetLineWidth());
1287  signalEst->SetFillStyle (backgroundEst->GetFillStyle());
1288  signalEst->SetFillColor (backgroundEst->GetFillColor());
1289  signalEst->SetTitle("Signal");
1290  signalEst->Add(backgroundEst,-1);
1291  bgCon->Add(signalEst);
1292 
1293  TH2* measured = GetH2(dataCon, "etaVsIPz");
1294  TH2* beta = static_cast<TH2*>(backgroundEst->Clone("beta"));
1295  beta->SetTitle("#beta");
1296  beta->SetDirectory(0);
1297  beta->Divide(measured);
1298  bgCon->Add(beta);
1299 
1300  TH2* bg1MBeta = static_cast<TH2*>(beta->Clone("oneMinusBeta"));
1301  bg1MBeta->SetTitle("1-#beta");
1302  bg1MBeta->SetDirectory(0);
1303  bg1MBeta->Reset();
1304  bgCon->Add(bg1MBeta);
1305  for (Int_t ix = 1; ix <= bg1MBeta->GetNbinsX(); ix++) {
1306  for (Int_t iy = 1; iy <= bg1MBeta->GetNbinsY(); iy++) {
1307  bg1MBeta->SetBinContent(ix,iy,1);
1308  bg1MBeta->SetBinError (ix,iy,0);
1309  }
1310  }
1311  bg1MBeta->Add(beta,-1);
1312 
1313  // if (lDebug > 2) bgCon->ls();
1314  return bgCon;
1315 }
1316 //____________________________________________________________________
1317 AliTrackletdNdetaTask::Container*
1318 AliTrackletdNdetaTask::CentBin::MasterFinalize(Container* parent,
1319  TH1* ipz,
1320  Double_t deltaCut,
1321  Double_t deltaTail)
1322 {
1323  Double_t nEvents = fIPz->GetEntries();
1324  Container* result = CreateContainer(*parent);
1325  Printf("Bin %5.1f-%5.1f%%: %10d events", fMin, fMax, Int_t(nEvents));
1326 
1327  // Copy ipZ histogram and scale by number of events
1328  TH1* ipZ = static_cast<TH1*>(CloneAndAdd(result, fIPz));
1329  ipZ->Scale(1./nEvents);
1330 
1331  Container* dataRes = 0;
1332  TIter next(fHistoSets);
1333  HistoSet* set = 0;
1334  while ((set = static_cast<HistoSet*>(next()))) {
1335  if (set == fDataSet) {
1336  dataRes = set->MasterFinalize(result, fIPz, deltaCut, deltaTail);
1337  continue;
1338  }
1339  EstimateBackground(dataRes, set, result, deltaCut, deltaTail);
1340  }
1341 
1342  if (lDebug > 5) {
1343  result->ls();
1344  }
1345  return result;
1346 }
1347 //____________________________________________________________________
1348 Bool_t AliTrackletdNdetaTask::CentBin::FillSet(UShort_t set,
1349  Bool_t isSignal,
1350  Double_t ipZ,
1351  Double_t eta,
1352  Double_t dPhi,
1353  Double_t dPhiS,
1354  Double_t dTheta,
1355  Double_t dThetaX,
1356  Double_t delta,
1357  Double_t weight)
1358 {
1359  HistoSet* hset = 0;
1360  if (set == kNormal) hset = fDataSet;
1361  else if (set == kInjection) hset = fInjSet;
1362  else if (set == kRotation) hset = fRotSet;
1363 
1364  if (!hset) return false;
1365 
1366  hset->Fill(isSignal ? ipZ : -1000,eta,dPhi,dPhiS,dTheta,dThetaX,delta,weight);
1367 
1368  return true;
1369 }
1370 
1371 //====================================================================
1372 // Task member functions
1373 //____________________________________________________________________
1374 AliTrackletdNdetaTask::AliTrackletdNdetaTask(const char* name)
1375  : AliAnalysisTaskSE(name),
1376  AliTrackletdNdetaUtils(),
1377  fContainer(0),
1378  fCentBins(0),
1379  fCentMethod("V0M"),
1380  fCentAxis(1,0,100),
1381  fIPzAxis(20,-10,10),
1382  fEtaAxis(80,-2,2),
1383  fPhiAxis(200,0,TMath::TwoPi()),
1384  fCent(0),
1385  fIPz(0),
1386  fEtaPhi(0),
1387  fUsedClusters0(0),
1388  fUsedClusters1(0),
1389  fAllClusters0(0),
1390  fAllClusters1(0),
1391  fStatus(0),
1392  fRecMode(kNormal|kInjection),
1393  fScaleDTheta(false),
1394  fDPhiShift(0.0045),
1395  fShiftedDPhiCut(-1),
1396  fScaledDThetaCut(-1),
1397  fDeltaCut(1.5),
1398  fMaxDelta(25),
1399  fTailDelta(5),
1400  fDThetaWindow(0.025),
1401  fDPhiWindow(0.06),
1402  fPhiOverlapCut(0.005),
1403  fZEtaOverlapCut(0.05),
1404  fPhiRotation(TMath::Pi())
1405 {
1406  FixAxis(fCentAxis, "Centrality [%]");
1407  FixAxis(fIPzAxis, "IP_{#it{z}} [cm]");
1408  FixAxis(fEtaAxis, "#eta");
1409  FixAxis(fPhiAxis, "#phi");
1410 
1411  // Set branches we want to look for
1412  fBranchNames =
1413  "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
1414  "SPDVertex.,PrimaryVertex.";
1415 
1416  DefineOutput(1, Container::Class());
1417  DefineOutput(2, Container::Class());
1418 }
1419 //____________________________________________________________________
1420 void AliTrackletdNdetaTask::UserCreateOutputObjects()
1421 {
1422  WorkerInit();
1423  lDebug = DebugLevel();
1424  PostData(1, fContainer);
1425 }
1426 //____________________________________________________________________
1427 void AliTrackletdNdetaTask::UserExec(Option_t*)
1428 {
1429  Double_t cent = -1;
1430  const AliVVertex* ip = 0;
1431  AliMultiplicity* mult = 0;
1432  TTree* clusters = 0;
1433  if (!CheckEvent(cent, ip, mult, clusters)) {
1434  AliWarningF("Event didn't pass %f, %p, %p, %p",
1435  cent, ip, mult, clusters);
1436  return;
1437  }
1438 
1439  ProcessEvent(cent, ip, mult, clusters);
1440 
1441  PostData(1,fContainer);
1442 
1443  fStatus->Fill(kCompleted);
1444 }
1445 //____________________________________________________________________
1446 void AliTrackletdNdetaTask::Terminate(Option_t*)
1447 {
1448  lDebug = DebugLevel();
1449 
1450  Container* results = new Container;
1451  results->SetName(Form("%sResults",GetName()));
1452  results->SetOwner();
1453 
1454  Print("");
1455  fContainer = static_cast<Container*>(GetOutputData(1));
1456  if (!fContainer) {
1457  AliWarning("No sum container found!");
1458  return;
1459  }
1460 
1461  InitCentBins(fContainer);
1462 
1463  MasterFinalize(results);
1464 
1465  PostData(2, results);
1466 }
1467 //____________________________________________________________________
1468 Bool_t AliTrackletdNdetaTask::Connect(const char* sumFile,
1469  const char* resFile)
1470 {
1471  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1472  if (!mgr) {
1473  AliError("No analysis manager to connect to.");
1474  return false;
1475  }
1476 
1477  // Add to the manager
1478  mgr->AddTask(this);
1479 
1480  // Create and connect output containers
1481  TString sumOut;
1482  TString resOut;
1483  if (sumFile && sumFile[0] != '\0') sumOut = sumFile;
1484  if (resFile && resFile[0] != '\0') resOut = resFile;
1485  else if (sumFile && sumFile[0] != '\0') resOut = sumFile;
1486 
1487  // If the string is null or 'default' connect to standard output file
1488  if (sumOut.IsNull() || sumOut.EqualTo("default", TString::kIgnoreCase))
1489  sumOut = AliAnalysisManager::GetCommonFileName();
1490  // If the string is null or 'default' connect to standard output file
1491  if (resOut.IsNull() || resOut.EqualTo("default", TString::kIgnoreCase))
1492  resOut = AliAnalysisManager::GetCommonFileName();
1493 
1494  // Always connect input
1495  mgr->ConnectInput(this, 0, mgr->GetCommonInputContainer());
1496 
1497  // Connect sum list unless the output 'none' is specified
1498  if (!sumOut.EqualTo("none", TString::kIgnoreCase)) {
1499  TString sumName(Form("%sSums", GetName()));
1500  AliAnalysisDataContainer* sumCon =
1501  mgr->CreateContainer(sumName, TList::Class(),
1502  AliAnalysisManager::kOutputContainer, sumOut);
1503  mgr->ConnectOutput(this, 1, sumCon);
1504  }
1505 
1506  // Connect the result list unless the output 'none' is specified
1507  if (!resOut.EqualTo("none", TString::kIgnoreCase)) {
1508  TString resName(Form("%sResults", GetName()));
1509  AliAnalysisDataContainer* resCon =
1510  mgr->CreateContainer(resName, TList::Class(),
1511  AliAnalysisManager::kParamContainer, resOut);
1512  mgr->ConnectOutput(this, 2, resCon);
1513  }
1514  return true;
1515 }
1516 //____________________________________________________________________
1517 void AliTrackletdNdetaTask::Print(Option_t*) const
1518 {
1519  Double_t shiftedDPhiCut = fShiftedDPhiCut;
1520  if (shiftedDPhiCut < 0)
1521  shiftedDPhiCut = TMath::Sqrt(fDeltaCut)*fDPhiWindow;
1522 
1523  Printf("%s: %s", ClassName(), GetName());
1524  Printf(" %22s: 0x%x", "Reconstruction mode", fRecMode);
1525  Printf(" %22s: %d", "Scale by sin^2(theta)", fScaleDTheta);
1526  Printf(" %22s: %f", "Delta phi shift", fDPhiShift);
1527  Printf(" %22s: %f", "Shifted Delta phi cut", shiftedDPhiCut);
1528  Printf(" %22s: %f", "Scaled Delta Theta cut", fScaledDThetaCut);
1529  Printf(" %22s: %f", "Delta cut", fDeltaCut);
1530  Printf(" %22s: %f", "max Delta", fMaxDelta);
1531  Printf(" %22s: %f", "tail Delta", fTailDelta);
1532  Printf(" %22s: %f", "Delta theta window", fDThetaWindow);
1533  Printf(" %22s: %f", "Delta phi window", fDPhiWindow);
1534  Printf(" %22s: %f", "phi overlap cut", fPhiOverlapCut);
1535  Printf(" %22s: %f", "z-eta overlap cut", fZEtaOverlapCut);
1536  Printf(" %22s: %f", "phi rotation", fPhiRotation);
1537  PrintAxis(fEtaAxis);
1538  PrintAxis(fPhiAxis);
1539  PrintAxis(fIPzAxis,1,"IPz");
1540  PrintAxis(fCentAxis,0);
1541 }
1542 //____________________________________________________________________
1543 void AliTrackletdNdetaTask::InitCentBins(Container* existing)
1544 {
1545  if (fCentBins) return;
1546  fCentBins = new Container;
1547  fCentBins->SetName("centralityBins");
1548  fCentBins->SetOwner();
1549 
1550  Double_t maxdPhi = fDPhiWindow *TMath::Sqrt(fMaxDelta);
1551  Double_t maxdTheta = fDThetaWindow*TMath::Sqrt(fMaxDelta);
1552  TAxis deltaAxis (Int_t(5*fMaxDelta),0, fMaxDelta);
1553  TAxis dThetaAxis(100, -maxdTheta,+maxdTheta);
1554  TAxis dPhiAxis (100, -maxdPhi, +maxdPhi);
1555  FixAxis(deltaAxis,
1556  Form("#Delta=[(#Delta#phi-#delta#phi)/#sigma_{#phi}]^{2}+"
1557  "[#Delta#theta%s/#sigma_{#theta}]^{2}",
1558  fScaleDTheta ? "sin^{-2}(#theta)" : ""));
1559  FixAxis(dThetaAxis, "#Delta#theta");
1560  FixAxis(dPhiAxis, "#Delta#phi");
1561 
1562  // Add min-bias bin
1563  CentBin* bin = MakeCentBin(0, 100, fRecMode);
1564  if (!existing)
1565  bin->WorkerInit(*fContainer,fEtaAxis,fIPzAxis,
1566  deltaAxis, dThetaAxis, dPhiAxis);
1567  else
1568  bin->FinalizeInit(existing);
1569  fCentBins->AddAt(bin, 0);
1570 
1571  // Add other bins
1572  Int_t nCentBins = fCentAxis.GetNbins();
1573  for (Int_t i = 1; i <= nCentBins; i++) {
1574  Float_t c1 = fCentAxis.GetBinLowEdge(i);
1575  Float_t c2 = fCentAxis.GetBinUpEdge(i);
1576  bin = MakeCentBin(c1, c2, fRecMode);
1577  if (!existing)
1578  bin->WorkerInit(*fContainer,fEtaAxis,fIPzAxis,
1579  deltaAxis, dThetaAxis, dPhiAxis);
1580  else
1581  bin->FinalizeInit(existing);
1582  fCentBins->AddAt(bin, i);
1583  }
1584 }
1585 //____________________________________________________________________
1586 Bool_t AliTrackletdNdetaTask::InitCDB()
1587 {
1588  Printf("Now initialising CDB");
1589  AliAnalysisManager* anaMgr = AliAnalysisManager::GetAnalysisManager();
1590  if (!anaMgr) {
1591  AliError("No manager defined!");
1592  return false;
1593  }
1594  // Check if we have the CDB connect task, and if so, do nothing as
1595  // we rely on that task set up things properly.
1596  const char* cdbNames[] = {"CDBconnect", "cdb", 0 };
1597  const char** ptr = cdbNames;
1598  while (*ptr) {
1599  AliAnalysisTask* cdbConnect = anaMgr->GetTask(*ptr);
1600  if (cdbConnect && cdbConnect->IsA()->InheritsFrom("AliTaskCDBconnect")) {
1601  AliInfoF("CDB-connect task (%s: %s) present, do nothing",
1602  cdbConnect->ClassName(), *ptr);
1603  return true;
1604  }
1605  ptr++;
1606  }
1607  // Otherwise, we need to do stuff ourselves
1608  Printf("Get the CDB manager");
1609  AliCDBManager* cdbMgr = AliCDBManager::Instance();
1610  if (!cdbMgr) {
1611  AliError("Failed to get instance of CDB manager");
1612  return false;
1613  }
1614  Int_t refRun = GetCDBReferenceRun();
1615  TString refUrl = GetCDBReferenceURL();
1616  AliWarningF("Using reference CDB storage \"%s\" and run \"%d\"",
1617  refUrl.Data(), refRun);
1618  // Set a very particular default storage. Perhaps we can do this
1619  // with specific storages instead! Depends on whether the
1620  // reconstruction also uses CDB - probably does - in which case we
1621  // need specific storages for that too.
1622  cdbMgr->SetDefaultStorage(refUrl);
1623  // Now load our geometry - from a LHC10h run
1624  Printf("Get Geometry entry");
1625  AliCDBEntry* cdbEnt = cdbMgr->Get("GRP/Geometry/Data", refRun);
1626  if (!cdbEnt) {
1627  AliErrorF("No geometry found from %d", refRun);
1628  return false;
1629  }
1630  // Initialize the geometry manager
1631  Printf("Set Geometry");
1632  AliGeomManager::SetGeometry(static_cast<TGeoManager*>(cdbEnt->GetObject()));
1633  // Now perform mis-alignment - again based on an LHC10h run!
1634  Printf("Misalign geometry");
1635  if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",refRun,-1,-1)) {
1636  AliErrorF("Failed to misalign geometry from %d", refRun);
1637  return false;
1638  }
1639  return true;
1640 }
1641 //____________________________________________________________________
1642 Bool_t AliTrackletdNdetaTask::InitGeometry()
1643 {
1644  if (fRecMode == 0) return true;
1645 
1646  if (!AliGeomManager::GetGeometry()) {
1647  AliError("No geometry loaded, needed for reconstruction");
1648  AliError("Add the AliTaskCDBconnect to the train");
1649  return false;
1650  }
1651  return true;
1652 }
1653 //____________________________________________________________________
1654 void AliTrackletdNdetaTask::WorkerInit()
1655 {
1656  if (!InitCDB()) {
1657  // Argh! - make us a a zombie
1658  SetZombie(true);
1659  return;
1660  }
1661 
1662  if (fShiftedDPhiCut < 0)
1663  fShiftedDPhiCut = TMath::Sqrt(fDeltaCut)*fDPhiWindow;
1664 
1665  fContainer = new Container;
1666  fContainer->SetName(Form("%sSums",GetName()));
1667  fContainer->SetOwner();
1668 
1669  fIPz = Make1D(*fContainer, "ipz", "", kMagenta+2, 20, fIPzAxis);
1670  fCent = Make1D(*fContainer, "cent", "", kMagenta+2, 20, fCentAxis);
1671  fEtaPhi = Make2D(*fContainer, "etaPhi","",kMagenta+2, 20, fEtaAxis,fPhiAxis);
1672  fIPz ->SetFillStyle(0);
1673  fCent ->SetFillStyle(0);
1674 
1675  TAxis zAxis(64,-16,16); FixAxis(zAxis, "#it{z} [cm]");
1676  TAxis phiAxis(80, 0,TMath::TwoPi()); FixAxis(phiAxis,"#phi");
1677  TAxis ph2Axis(160,0,TMath::TwoPi()); FixAxis(ph2Axis,"#phi");
1678 
1679  fUsedClusters0 = Make2D(*fContainer, "usedClusters0",
1680  "Used clusters on layer 0",
1681  kBlack, 1, zAxis, phiAxis);
1682  fUsedClusters1 = Make2D(*fContainer, "usedClusters1",
1683  "Used clusters on layer 1",
1684  kBlack, 1, zAxis, ph2Axis);
1685  fAllClusters0 = Make2D(*fContainer, "allClusters0",
1686  "All clusters on layer 0",
1687  kBlack, 1, zAxis, phiAxis);
1688  fAllClusters1 = Make2D(*fContainer, "allClusters1",
1689  "All clusters on layer 1",
1690  kBlack, 1, zAxis, ph2Axis);
1691 
1692  fStatus = new TH1F("status", "Status of task",
1693  kCompleted, .5, kCompleted+.5);
1694  fStatus->SetMarkerSize(2);
1695  fStatus->SetMarkerColor(kMagenta+2);
1696  fStatus->SetLineColor(kMagenta+2);
1697  fStatus->SetFillColor(kMagenta+2);
1698  fStatus->SetFillStyle(1001);
1699  fStatus->SetBarOffset(0.1);
1700  fStatus->SetBarWidth(0.4);
1701  fStatus->SetDirectory(0);
1702  fStatus->SetStats(0);
1703  fStatus->SetXTitle("Event have");
1704  fStatus->SetYTitle("# Events");
1705  fStatus->GetXaxis()->SetBinLabel(kAll, "Been seen");
1706  fStatus->GetXaxis()->SetBinLabel(kEvent, "Event data");
1707  fStatus->GetXaxis()->SetBinLabel(kMultiplicity, "Multiplicity");
1708  fStatus->GetXaxis()->SetBinLabel(kClusters, "Clusters");
1709  fStatus->GetXaxis()->SetBinLabel(kTrigger, "Trigger");
1710  fStatus->GetXaxis()->SetBinLabel(kIP, "IP");
1711  fStatus->GetXaxis()->SetBinLabel(kCentrality, "Centrality");
1712  fStatus->GetXaxis()->SetBinLabel(kReconstructed, "Been reconstructed");
1713  fStatus->GetXaxis()->SetBinLabel(kInjected, "Done Injection");
1714  fStatus->GetXaxis()->SetBinLabel(kRotated, "Been rotated");
1715  fStatus->GetXaxis()->SetBinLabel(kCompleted, "Completed");
1716  fContainer->Add(fStatus);
1717 
1718  typedef TParameter<double> DP;
1719  typedef TParameter<bool> BP;
1720  typedef TParameter<int> IP;
1721  Container* params = new Container;
1722  params->SetName("parameters");
1723  params->SetOwner();
1724  fContainer->Add(params);
1725  params->Add(new IP("RecMode", fRecMode, 'f'));
1726  params->Add(new BP("ScaleDTheta", fScaleDTheta, 'f'));
1727  params->Add(new DP("DPhiShift", fDPhiShift, 'f'));
1728  params->Add(new DP("ShiftedDPhiCut", fShiftedDPhiCut, 'f'));
1729  params->Add(new DP("ScaledDThetaCut",fScaledDThetaCut,'f'));
1730  params->Add(new DP("DeltaCut", fDeltaCut, 'f'));
1731  params->Add(new DP("MaxDelta", fMaxDelta, 'f'));
1732  params->Add(new DP("TailDelta", fTailDelta, 'f'));
1733  params->Add(new DP("DThetaWindow", fDThetaWindow, 'f'));
1734  params->Add(new DP("DPhiWindow", fDPhiWindow, 'f'));
1735  params->Add(new DP("PhiOverlapCut", fPhiOverlapCut, 'f'));
1736  params->Add(new DP("ZEtaOverlapCut" ,fZEtaOverlapCut, 'f'));
1737  params->Add(new DP("PhiRotation", fPhiRotation, 'f'));
1738 
1739  // Fix some fill styles
1740 
1741  InitCentBins(0);
1742 }
1743 //____________________________________________________________________
1744 void AliTrackletdNdetaTask::WorkerFinalize()
1745 {
1746  Printf("This worker had %f events within cuts",
1747  fIPz->GetEntries());
1748 }
1749 //____________________________________________________________________
1750 void AliTrackletdNdetaTask::MasterFinalize(Container* results)
1751 {
1752  // Make copies of histograms and store
1753  fIPz = static_cast<TH1*>(CloneAndAdd(results,GetH1(fContainer,"ipz")));
1754  fStatus = static_cast<TH1*>(CloneAndAdd(results,GetH1(fContainer,"status")));
1755  fCent = static_cast<TH1*>(CloneAndAdd(results,GetH1(fContainer,"cent")));
1756 
1757  Double_t nEvents = fIPz->GetEntries();
1758  Printf("Event summary:");
1759  for (Int_t i = 1; i <= fStatus->GetNbinsX(); i++)
1760  Printf(" %10d %s", Int_t(fStatus->GetBinContent(i)),
1761  fStatus->GetXaxis()->GetBinLabel(i));
1762 
1763  fIPz ->Scale(1./nEvents);
1764  fStatus->Scale(1./fStatus->GetBinContent(1));
1765  fCent ->Scale(1./fCent->GetEntries());
1766 
1767  fUsedClusters0 = static_cast<TH2*>(CloneAndAdd(results,
1768  GetH2(fContainer,
1769  "usedClusters0")));
1770  fUsedClusters1 = static_cast<TH2*>(CloneAndAdd(results,
1771  GetH2(fContainer,
1772  "usedClusters1")));
1773  fAllClusters0 = static_cast<TH2*>(CloneAndAdd(results,
1774  GetH2(fContainer,
1775  "allClusters0")));
1776  fAllClusters1 = static_cast<TH2*>(CloneAndAdd(results,
1777  GetH2(fContainer,
1778  "allClusters1")));
1779  fUsedClusters0->Scale(1/nEvents);
1780  fUsedClusters1->Scale(1/nEvents);
1781  fAllClusters0 ->Scale(1/nEvents);
1782  fAllClusters1 ->Scale(1/nEvents);
1783 
1784  TIter next(fCentBins);
1785  CentBin* bin = 0;
1786  while ((bin = static_cast<CentBin*>(next()))) {
1787  bin->MasterFinalize(results, 0, fDeltaCut, fTailDelta);
1788  }
1789 }
1790 //____________________________________________________________________
1791 TTree* AliTrackletdNdetaTask::FindClusters(Bool_t needed)
1792 {
1793  if (!needed) return 0;
1794  AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
1795  AliVEventHandler* inh = mgr->GetInputEventHandler();
1796  if (!inh->IsA()->InheritsFrom(AliESDInputHandlerRP::Class())) {
1797  AliErrorF("Not the right kind of input handler: %s",
1798  inh->ClassName());
1799  return 0;
1800  }
1801  AliESDInputHandlerRP* rph = static_cast<AliESDInputHandlerRP*>(inh);
1802  TTree* tree = rph->GetTreeR("ITS");
1803  if (!tree) {
1804  AliError("Tree of clusters (rec.points) not found");
1805  return 0;
1806  }
1807  return tree;
1808 }
1809 //____________________________________________________________________
1810 Double_t AliTrackletdNdetaTask::FindCentrality(AliVEvent* event)
1811 {
1812  if (fCentMethod.EqualTo("MB", TString::kIgnoreCase)) return 0;
1813  AliMultSelection* cent =
1814  static_cast<AliMultSelection*>(event->FindListObject("MultSelection"));
1815  if (!cent) {
1816  AliWarning("No centrality in event");
1817  return -1;
1818  }
1819  const Double_t safety = 1e-3;
1820  Double_t centPer = cent->GetMultiplicityPercentile(fCentMethod);
1821  if (centPer < -safety) return -2;
1822  if (centPer < +safety) centPer = safety;
1823  else if (centPer > 100-safety) centPer = 100-safety;
1824 
1825  if (centPer < fCentAxis.GetXmin() || centPer > fCentAxis.GetXmax()) {
1826  AliWarningF("Centrality = %f out of range [%f,%f]",
1827  centPer, fCentAxis.GetXmin(), fCentAxis.GetXmax());
1828  return -3;
1829  }
1830  return centPer;
1831 }
1832 //____________________________________________________________________
1833 const AliVVertex* AliTrackletdNdetaTask::FindIP(AliVEvent* event,
1834  Double_t maxDispersion,
1835  Double_t maxZError)
1836 {
1837  const AliVVertex* ip = event->GetPrimaryVertex();
1838  if (ip->GetNContributors() <= 0) {
1839  AliWarning("Not enough contributors for IP");
1840  return 0;
1841  }
1842  // If this is from the Z vertexer, do some checks
1843  if (ip->IsFromVertexerZ()) {
1844  // Get covariance matrix
1845  Double_t covar[6];
1846  ip->GetCovarianceMatrix(covar);
1847  Double_t sigmaZ = TMath::Sqrt(covar[5]);
1848  if (sigmaZ >= maxZError) {
1849  AliWarningF("IPz resolution = %f >= %f", sigmaZ, maxZError);
1850  return 0;
1851  }
1852 
1853  // If this IP doesn not derive from AliVertex, don't check dispersion.
1854  if (ip->IsA()->InheritsFrom(AliVertex::Class())) {
1855  const AliVertex* ipv = static_cast<const AliVertex*>(ip);
1856  // Dispersion is the parameter used by the vertexer for finding the IP.
1857  if (ipv->GetDispersion() >= maxDispersion) {
1858  AliWarningF("IP dispersion = %f >= %f",
1859  ipv->GetDispersion(), maxDispersion);
1860  return 0;
1861  }
1862  }
1863  }
1864 
1865  // If we get here, we either have a full 3D vertex or track
1866  // vertex, and we should check if it is in range
1867  if (ip->GetZ() < fIPzAxis.GetXmin() || ip->GetZ() > fIPzAxis.GetXmax()) {
1868  AliWarningF("IPz = %fcm out of range [%f,%f]cm",
1869  ip->GetZ(), fIPzAxis.GetXmin(), fIPzAxis.GetXmax());
1870  return 0;
1871  }
1872  // Good vertex, return it
1873  return ip;
1874 }
1875 //____________________________________________________________________
1876 Bool_t AliTrackletdNdetaTask::FindTrigger()
1877 {
1878  UInt_t evBits = fInputHandler->IsEventSelected();
1879  Bool_t trgOK = evBits & fOfflineTriggerMask;
1880 
1881  return trgOK;
1882 }
1883 //____________________________________________________________________
1884 Bool_t AliTrackletdNdetaTask::ProcessTracklets(TList& toRun,
1885  AliITSMultRecBg* reco,
1886  Double_t cent,
1887  Double_t ipZ,
1888  AliMultiplicity* mult)
1889 {
1890  // This can probably be optimized by using the reco object
1891  // directly. For now, we leave as is and accept the extra
1892  // processing time it implies.
1893  if (!reco) FillClusters(mult, ipZ);
1894  else if (FindMode(reco) == kNormal) FillClusters(reco);
1895 
1896  // Now loop over all tracklets. Since this method is only called
1897  // once per "reconstruction" pass, we save a lot of computing time
1898  // since we loop over the tracklets of each reconstruction exactly
1899  // once. It also means that we implement the calculations and
1900  // cuts in exacly one place (except caveat on filling clusters above).
1901  Int_t nTracklets = mult->GetNumberOfTracklets();
1902  for (Int_t trackletNumber = 0; trackletNumber < nTracklets;
1903  trackletNumber++) {
1904  Double_t weight = TrackletWeight(trackletNumber, mult);
1905  Double_t theta = mult->GetTheta(trackletNumber);
1906  Double_t phi = mult->GetPhi(trackletNumber);
1907  Double_t dTheta = mult->GetDeltaTheta(trackletNumber);
1908  Double_t dPhi = mult->GetDeltaPhi(trackletNumber);
1909  Double_t delta = mult->CalcDist(trackletNumber);
1910  Double_t eta = -TMath::Log(TMath::Tan(theta/2));
1911  Double_t dThetaX = dTheta;
1912  Double_t dPhiS = dPhi - TMath::Sign(fDPhiShift,dPhi); //+ if dPhi<0
1913  if (fScaleDTheta) dThetaX /= TMath::Power(TMath::Sin(theta),2);
1914 
1915  // Optional cut on scaled DTheta
1916  if (fScaledDThetaCut >= 0 &&
1917  TMath::Abs(dThetaX) > fScaledDThetaCut) continue;
1918  // If not within the bins we're looking at, continue
1919  if (eta < fEtaAxis.GetXmin() || eta > fEtaAxis.GetXmax()) continue;
1920  // if (phi < fPhiAxis.GetXmin() || phi > fPhiAxis.GetXmax()) continue;
1921 
1922  // Assume a signal
1923  Bool_t isSignal = true;
1924  if (delta > fDeltaCut) isSignal = false;
1925  if (dPhiS >= fShiftedDPhiCut) isSignal = false;
1926 
1927  // Now fill all found centrality bins
1928  FillBins(toRun, reco, mult, trackletNumber, isSignal, ipZ,
1929  eta, dPhi, dPhiS, dTheta, dThetaX, delta, weight);
1930 
1931  // If this is a signal, fill into eta,phi map
1932  if (isSignal) fEtaPhi->Fill(eta,phi);
1933  }
1934  return true; // nAcc > 0;
1935 }
1936 //____________________________________________________________________
1937 UShort_t AliTrackletdNdetaTask::FindMode(AliITSMultRecBg* reco) const
1938 {
1939  if (!reco) return kNormal;
1940  switch (reco->GetRecType()) {
1941  case AliITSMultRecBg::kData: return kNormal;
1942  case AliITSMultRecBg::kBgInj: return kInjection;
1943  case AliITSMultRecBg::kBgRot: return kRotation;
1944  default:
1945  AliWarningF("Unsupported reconstruction mode: %d", reco->GetRecType());
1946  break;
1947  }
1948  return 999;
1949 }
1950 //____________________________________________________________________
1951 void AliTrackletdNdetaTask::FillBins(TList& toRun,
1952  AliITSMultRecBg* reco,
1953  AliMultiplicity* mult,
1954  Int_t trackletNumber,
1955  Bool_t isSignal,
1956  Double_t ipz,
1957  Double_t eta,
1958  Double_t dPhi,
1959  Double_t dPhiS,
1960  Double_t dTheta,
1961  Double_t dThetaX,
1962  Double_t delta,
1963  Double_t weight) const
1964 {
1965  UShort_t mode = FindMode(reco);
1966  if (mode >= 999) return;
1967 
1968  TIter nextB(&toRun);
1969  CentBin* bin = 0;
1970  while ((bin = static_cast<CentBin*>(nextB())))
1971  bin->FillSet(mode, isSignal, ipz, eta,
1972  dPhi, dPhiS, dTheta, dThetaX, delta, weight);
1973 }
1974 //____________________________________________________________________
1975 Bool_t AliTrackletdNdetaTask::Reconstruct(TList& toRun,
1976  UShort_t mode,
1977  TTree* clusters,
1978  Double_t cent,
1979  const AliVVertex* ip)
1980 {
1981  UShort_t firstMode = mode;
1982  AliITSMultRecBg reco;
1983  reco.SetCreateClustersCopy (true);
1984  reco.SetScaleDThetaBySin2T (fScaleDTheta);
1985  reco.SetNStdDev (fMaxDelta);
1986  reco.SetPhiWindow (fDPhiWindow);
1987  reco.SetThetaWindow (fDThetaWindow);
1988  reco.SetPhiShift (fDPhiShift);
1989  reco.SetRemoveClustersFromOverlaps(fPhiOverlapCut > 0 ||
1990  fZEtaOverlapCut > 0);
1991  reco.SetPhiOverlapCut (fPhiOverlapCut);
1992  reco.SetZetaOverlapCut (fZEtaOverlapCut);
1993  reco.SetHistOn (false);
1994  switch (mode) {
1995  case kNormal:
1996  case kInjection:
1997  // If we're to do injection, first run normally, since the
1998  // injection procedure relies on this, and there's no reason to
1999  // do it twice
2000  firstMode = kNormal;
2001  reco.SetRecType(AliITSMultRecBg::kData);
2002  break;
2003  case kRotation:
2004  reco.SetRecType(AliITSMultRecBg::kBgRot);
2005  reco.SetPhiRotationAngle(fPhiRotation);
2006  break;
2007  }
2008 
2009  // Now run the reconstruction
2010  Float_t ipv[3];
2011  ipv[0] = ip->GetX();
2012  ipv[1] = ip->GetY();
2013  ipv[2] = ip->GetZ();
2014  reco.Run(clusters, ipv);
2015 
2016  // And fill results into centrality bins
2017  Bool_t ret = ProcessTracklets(toRun, &reco, cent, ip->GetZ(),
2018  reco.GetMultiplicity());
2019 
2020  // If we're not doing injection, return
2021  if (mode != kInjection) return ret;
2022 
2023  // If we do injection, run again, but in injection mode
2024  reco.SetRecType(AliITSMultRecBg::kBgInj);
2025  reco.Run(clusters, ipv);
2026 
2027  // And fill results into centrality bins
2028  ret = ProcessTracklets(toRun, &reco, cent, ip->GetZ(),
2029  reco.GetMultiplicity());
2030 
2031  return ret;
2032 }
2033 //____________________________________________________________________
2034 void AliTrackletdNdetaTask::FillClusters(AliITSMultRecBg* reco)
2035 {
2036  if (!reco) return;
2037 
2038  TH2* usedClusters[] = { fUsedClusters0, fUsedClusters1 };
2039  TH2* allClusters[] = { fAllClusters0, fAllClusters1 };
2040 
2041  // Loop over tracklets
2042  Int_t nTracklets = reco->GetNTracklets();
2043  for (Int_t trackletNo = nTracklets; trackletNo--; ) {
2044  // Get tracklet parameters
2045  Float_t* tracklet = reco->GetTracklet(trackletNo);
2046 
2047  Float_t dPhi = tracklet[AliITSMultReconstructor::kTrDPhi];
2048  Float_t dTheta = tracklet[AliITSMultReconstructor::kTrDTheta];
2049  Float_t theta = tracklet[AliITSMultReconstructor::kTrTheta];
2050  if (TMath::Abs(TMath::Abs(dPhi)-fDPhiShift) > fShiftedDPhiCut)
2051  continue;
2052 
2053  Float_t delta = reco->CalcDist(dPhi, dTheta, theta);
2054  if (delta > fDeltaCut) continue;
2055 
2056  Int_t clusterID[] = { Int_t(tracklet[AliITSMultReconstructor::kClID1]),
2057  Int_t(tracklet[AliITSMultReconstructor::kClID2]) };
2058  for (Int_t layer = 0; layer < 2; layer++) {
2059  // Get cluster information
2060  Float_t* cluster = reco->GetClusterOfLayer(layer, clusterID[layer]);
2061  // Fill used cluster z,phi
2062  usedClusters[layer]->Fill(cluster[AliITSMultReconstructor::kClZ],
2063  cluster[AliITSMultReconstructor::kClPh]);
2064  }
2065  }
2066  // Loop over all clusters and fill information
2067  for (Int_t layer = 0; layer < 2; layer++) {
2068  for (Int_t clusterNo = reco->GetNClustersLayer(layer); clusterNo--; ) {
2069  Float_t* cluster = reco->GetClusterOfLayer(layer, clusterNo);
2070  allClusters[layer]->Fill(cluster[AliITSMultReconstructor::kClZ],
2071  cluster[AliITSMultReconstructor::kClPh]);
2072  }
2073  }
2074 }
2075 //____________________________________________________________________
2076 void AliTrackletdNdetaTask::FillClusters(AliMultiplicity* mult, Double_t ipZ)
2077 {
2078  if (!mult) return;
2079  const double kRSPD2 = 3.9;
2080 
2081  // Loop over tracklets
2082  Int_t nTracklets = mult->GetNumberOfTracklets();
2083  for (Int_t trackletNo = nTracklets; trackletNo--; ) {
2084  Double_t phi = mult->GetPhi(trackletNo);
2085  Double_t z = kRSPD2 / TMath::Tan(mult->GetTheta(trackletNo)) + ipZ;
2086  Double_t dPhi = mult->GetDeltaPhi(trackletNo);
2087  if ((TMath::Abs(dPhi-fDPhiShift) <= fShiftedDPhiCut) &&
2088  (mult->CalcDist(trackletNo )<= fDeltaCut))
2089  fUsedClusters0->Fill(z, phi);
2090  fAllClusters0->Fill(z, phi);
2091  }
2092  // Loop over free clusters
2093  Int_t nClusters = mult->GetNumberOfSingleClusters();
2094  for (Int_t clusterNo = nClusters; clusterNo--; ) {
2095  Double_t phi = mult->GetPhiSingle(clusterNo);
2096  Double_t z = kRSPD2 / TMath::Tan(mult->GetThetaSingle(clusterNo)) + ipZ;
2097  fAllClusters0->Fill(z, phi);
2098  }
2099 }
2100 //____________________________________________________________________
2101 Bool_t AliTrackletdNdetaTask::CheckEvent(Double_t& cent,
2102  const AliVVertex*& ip,
2103  AliMultiplicity*& mult,
2104  TTree*& clusters)
2105 {
2106  // Count all events
2107  fStatus->Fill(kAll);
2108 
2109  // Check for event
2110  AliVEvent* event = InputEvent();
2111  if (!event) {
2112  AliWarning("No event");
2113  return false;
2114  }
2115  fStatus->Fill(kEvent);
2116 
2117  // Check for geometry
2118  if (!InitGeometry()) {
2119  // Argh! - make us a a zombie
2120  SetZombie(true);
2121  return false;
2122  }
2123 
2124  // Check magnetic field
2125  if (!TGeoGlobalMagField::Instance()->GetField() &&
2126  !event->InitMagneticField()) {
2127  AliWarning("Failed to initialize magnetic field");
2128  return false;
2129  }
2130 
2131  // Check the multiplicity
2132  AliVMultiplicity* vmult = event->GetMultiplicity();
2133  if (!vmult) return false;
2134  if (!vmult->IsA()->InheritsFrom(AliMultiplicity::Class())) {
2135  AliWarningF("Multiplicity is a %s object", vmult->ClassName());
2136  return false;
2137  }
2138  mult = static_cast<AliMultiplicity*>(vmult);
2139  fStatus->Fill(kMultiplicity);
2140 
2141  // Check that we have clusters - if so needed
2142  Bool_t needRP = (fRecMode != 0);
2143  clusters = FindClusters(needRP);
2144  if (needRP && !clusters) return false;
2145  if (needRP) fStatus->Fill(kClusters);
2146 
2147  // Check if event was triggered
2148  Bool_t trg = FindTrigger();
2149  if (!trg) return false;
2150  fStatus->Fill(kTrigger);
2151 
2152  // Check the interaction point
2153  ip = FindIP(event);
2154  if (!ip) return false;
2155  fStatus->Fill(kIP);
2156 
2157  // Check the centrality
2158  cent = FindCentrality(event);
2159  if (cent < 0) return false;
2160  fStatus->Fill(kCentrality);
2161 
2162  fIPz->Fill(ip->GetZ());
2163  fCent->Fill(cent);
2164  return true;
2165 }
2166 
2167 //____________________________________________________________________
2168 void AliTrackletdNdetaTask::ProcessEvent(Double_t cent,
2169  const AliVVertex* ip,
2170  AliMultiplicity* mult,
2171  TTree* clusters)
2172 {
2173  // Figure out which centrality bins to fill
2174  Int_t nAcc = 0;
2175  TIter next(fCentBins);
2176  CentBin* bin = 0;
2177  TList toRun;
2178  while ((bin = static_cast<CentBin*>(next()))) {
2179  if (!bin->Accept(cent)) continue; // Not in range for this bin
2180  toRun.Add(bin);
2181  bin->FillIPz(ip->GetZ());
2182  nAcc++;
2183  }
2184  // If we have no centrality bins to fill, we return immediately
2185  if (nAcc <= 0) return;
2186 
2187  if (!(fRecMode & (kNormal | kInjection))) {
2188  // In case we do not do normal reconstruction nor injection, we
2189  // should fill our data histograms from the read (ESD)
2190  // multiplicity.
2191  ProcessTracklets(toRun, 0, cent, ip->GetZ(), mult);
2192  }
2193  if ((fRecMode & kNormal) && !(fRecMode & kInjection)) {
2194  // Only do explicit normal reconstruction if we're not
2195  // injecting. In case of injection, we will also run the normal
2196  // mode with the same background estimator object.
2197  Reconstruct(toRun, kNormal, clusters, cent, ip);
2198  fStatus->Fill(kReconstructed);
2199  }
2200  if (fRecMode & kInjection) {
2201  Reconstruct(toRun, kInjection, clusters, cent, ip);
2202  fStatus->Fill(kReconstructed);
2203  fStatus->Fill(kInjected);
2204  }
2205  if (fRecMode & kRotation) {
2206  Reconstruct(toRun, kRotation, clusters, cent, ip);
2207  fStatus->Fill(kRotated);
2208  }
2209 }
2210 
2211 
2212 /*********************************************************************
2213  *
2214  * Code for processing simulations
2215  */
2216 #ifndef __CINT__
2217 #include <AliStack.h>
2218 #include <AliMCEvent.h>
2219 #include <AliGenEventHeader.h>
2220 #include <TParticle.h>
2221 #include <TParticlePDG.h>
2222 #include <TDatabasePDG.h>
2223 #include <TGraphAsymmErrors.h>
2224 #else
2225 class TParticlePDG; // Auto-load
2226 class TGraphAsymmErrors; // Auto-load
2227 class TParticle;
2228 class AliStack;
2229 #endif
2230 
2231 //====================================================================
2237 class AliTrackletdNdetaMCTask : public AliTrackletdNdetaTask
2238 {
2239 public:
2240  enum {
2241  kPrimaries = 0x8,
2242  kSecondaries = 0x10,
2243  kCombinatorics = 0x20,
2244  kUncorrCombinatorics = 0x40
2245  };
2246 
2247  //==================================================================
2248  struct CentBin : public AliTrackletdNdetaTask::CentBin
2249  {
2250  CentBin()
2251  : AliTrackletdNdetaTask::CentBin(),
2252  fPrimSet(0),
2253  fSecSet(0),
2254  fCombSet(0),
2255  fCombUSet(0),
2256  fEtaVsIPzMC(0),
2257  fEtaVsIPzMCSel(0),
2258  fPdgMC(0),
2259  fPrimaryPdg(0),
2260  fSecondaryPdg(0),
2261  fPrimaryParentPdg(0),
2262  fSecondaryParentPdg(0),
2263  fIPzVsMC(0),
2264  fIPzGen(0),
2265  fIPzSel(0)
2266  {}
2267  CentBin(Float_t c1, Float_t c2, UShort_t recFlags);
2279  virtual void WorkerInit(Container& parent,
2280  const TAxis& etaAxis,
2281  const TAxis& ipzAxis,
2282  const TAxis& deltaAxis,
2283  const TAxis& dThetaAxis,
2284  const TAxis& dPhiAxis);
2302  virtual Bool_t FillSet(UShort_t set,
2303  Bool_t isSignal,
2304  Double_t ipZ,
2305  Double_t eta,
2306  Double_t dPhi,
2307  Double_t dPhiS,
2308  Double_t dTheta,
2309  Double_t dThetaX,
2310  Double_t delta,
2311  Double_t weight);
2322  virtual void FillSpecie(Bool_t primary,
2323  Int_t pdg,
2324  Int_t parentPdg,
2325  Double_t delta,
2326  Double_t weight);
2333  virtual void FillMCIPz(Double_t ipZ, Double_t genIPz);
2344  virtual void FillPrimary(Double_t ipZ, Double_t genIPz,
2345  Double_t eta, Bool_t charged,
2346  Int_t pdg, Double_t weight);
2357  virtual TH2* MakeAlpha(const char* name,
2358  const char* title,
2359  Container* result,
2360  TH2* primaries) const;
2371  virtual TH2* MakeAlphaMask(TH2* a1, TH2* a2,
2372  Double_t min=0, Double_t max=2.5) const;
2383  virtual TH1* MakePdgProj(const char* name,
2384  const char* title,
2385  TH2* orig,
2386  Double_t cut=-1) const;
2397  virtual Container* MasterFinalize(Container* parent,
2398  TH1* ipz,
2399  Double_t deltaCut,
2400  Double_t deltaTail);
2401  protected:
2407  CentBin(const CentBin& o)
2408  : AliTrackletdNdetaTask::CentBin(o),
2409  fPrimSet(0),
2410  fSecSet(0),
2411  fCombSet(0),
2412  fCombUSet(0),
2413  fEtaVsIPzMC(0),
2414  fEtaVsIPzMCSel(0),
2415  fPdgMC(0),
2416  fPrimaryPdg(0),
2417  fSecondaryPdg(0),
2418  fPrimaryParentPdg(0),
2419  fSecondaryParentPdg(0),
2420  fIPzVsMC(0),
2421  fIPzGen(0),
2422  fIPzSel(0)
2423  {}
2430  CentBin& operator=(const CentBin&) { return *this; }
2431 
2432  HistoSet* fPrimSet;
2433  HistoSet* fSecSet;
2434  HistoSet* fCombSet;
2435  HistoSet* fCombUSet;
2436  TH2* fEtaVsIPzMC;
2437  TH2* fEtaVsIPzMCSel;
2438  TH2* fIPzVsMC;
2439  TH1* fIPzGen;
2440  TH1* fIPzSel;
2441  TH1* fPdgMC;
2442  TH2* fPrimaryPdg;
2443  TH2* fSecondaryPdg;
2444  TH2* fPrimaryParentPdg;
2445  TH2* fSecondaryParentPdg;
2446  ClassDef(CentBin,1); // Centrality bin
2447  };
2448  //==================================================================
2452  AliTrackletdNdetaMCTask()
2453  : AliTrackletdNdetaTask(),
2454  fMCStatus(0),
2455  fReweigh(false)
2456  {}
2460  AliTrackletdNdetaMCTask(const char* name, Bool_t reweigh=false)
2461  : AliTrackletdNdetaTask(name),
2462  fMCStatus(0),
2463  fReweigh(reweigh)
2464  {}
2473  virtual AliTrackletdNdetaTask::CentBin*
2474  MakeCentBin(Float_t c1, Float_t c2, UShort_t recFlags)
2475  {
2476  return new AliTrackletdNdetaMCTask::CentBin(c1, c2, recFlags);
2477  }
2478 protected:
2479 
2487  AliTrackletdNdetaMCTask(const AliTrackletdNdetaMCTask& o)
2488  : AliTrackletdNdetaTask(o),
2489  fMCStatus(0),
2490  fReweigh(false)
2491  {}
2497  AliTrackletdNdetaMCTask& operator=(const AliTrackletdNdetaMCTask& o)
2498  {
2499  return *this;
2500  }
2505  virtual void Print(Option_t* option) const
2506  {
2507  AliTrackletdNdetaTask::Print(option);
2508 
2509  Printf(" %22s: %s", "Reweigh", fReweigh ? "yes" : "no");
2510  }
2511  /* @} */
2516  enum {
2517  kAllMC = 1,
2518  kMCEvent,
2519  kStack,
2520  kGenHeader,
2521  kDataOK,
2522  kMCTruth,
2523  kReweighed,
2524  kMCCompleted
2525  };
2531  virtual const char* GetCDBReferenceURL() const
2532  {
2533  return "alien://Folder=/alice/simulation/2008/v4-15-Release/Residual";
2534  }
2539  virtual void WorkerInit();
2551  virtual Bool_t CheckEvent(Double_t& cent,
2552  const AliVVertex*& ip,
2553  AliMultiplicity*& mult,
2554  TTree*& clusters);
2560  virtual void MasterFinalize(Container* results);
2569  virtual Double_t TrackletWeight(Int_t trackletNumber,
2570  const AliMultiplicity* mult) const;
2578  void ReweighStack(Double_t cent) const;
2590  Double_t LookupWeight(Int_t trackNo,
2591  TParticle* particle,
2592  Double_t cent) const;
2603  virtual Double_t LookupWeight(TParticle* particle, Double_t cent) const
2604  {
2605  return 1;
2606  }
2615  virtual Bool_t FillPrimaries(Double_t cent,
2616  Double_t ipZ,
2617  Double_t genIPz,
2618  Bool_t dataOK);
2628  Bool_t HaveCommonParent(const Float_t* labels0,
2629  const Float_t* labels1) const;
2647  virtual void FillBins(TList& toRun,
2648  AliITSMultRecBg* reco,
2649  AliMultiplicity* mult,
2650  Int_t trackletNumber,
2651  Bool_t isSignal,
2652  Double_t ipz,
2653  Double_t eta,
2654  Double_t dPhi,
2655  Double_t dPhiS,
2656  Double_t dTheta,
2657  Double_t dThetaX,
2658  Double_t delta,
2659  Double_t weight) const;
2660  /* @} */
2672  static Int_t* PdgArray(Int_t& size);
2680  static Int_t* SortedPdgArray(Int_t& size);
2691  static Int_t FindPdg(Int_t pdg);
2698  static TAxis* MakePdgAxis();
2699  /* @} */
2700  TH1* fMCStatus; // Status of task
2701  Bool_t fReweigh;
2702  ClassDef(AliTrackletdNdetaMCTask,1); // MC class for dN/deta u/tracklets
2703 };
2704 
2705 //====================================================================
2706 AliTrackletdNdetaMCTask::CentBin::CentBin(Float_t c1, Float_t c2,
2707  UShort_t recFlags)
2708  : AliTrackletdNdetaTask::CentBin(c1,c2,recFlags, +2, 4),
2709  fPrimSet(0),
2710  fSecSet(0),
2711  fCombSet(0),
2712  fCombUSet(0),
2713  fEtaVsIPzMC(0),
2714  fEtaVsIPzMCSel(0),
2715  fPdgMC(0),
2716  fPrimaryPdg(0),
2717  fSecondaryPdg(0),
2718  fPrimaryParentPdg(0),
2719  fSecondaryParentPdg(0),
2720  fIPzVsMC(0),
2721  fIPzGen(0),
2722  fIPzSel(0)
2723 {
2724  fPrimSet = new HistoSet("primaries", kCyan+2, 34);
2725  fSecSet = new HistoSet("secondaries", kMagenta+2, 28);
2726  fCombSet = new HistoSet("combinatorics", kYellow+2, 30);
2727  fCombUSet = new HistoSet("uncorrelatedCombinatorics", kPink+2, 27);
2728  fHistoSets->Add(fPrimSet);
2729  fHistoSets->Add(fSecSet);
2730  fHistoSets->Add(fCombSet);
2731  fHistoSets->Add(fCombUSet);
2732 }
2733 //____________________________________________________________________
2734 void AliTrackletdNdetaMCTask::CentBin::WorkerInit(Container& parent,
2735  const TAxis& etaAxis,
2736  const TAxis& ipzAxis,
2737  const TAxis& deltaAxis,
2738  const TAxis& dThetaAxis,
2739  const TAxis& dPhiAxis)
2740 {
2741  AliTrackletdNdetaTask::CentBin::WorkerInit(parent, etaAxis,
2742  ipzAxis, deltaAxis,
2743  dThetaAxis, dPhiAxis);
2744  Container* c = fContainer;
2745  TAxis& pdg = *MakePdgAxis();
2746  fEtaVsIPzMC = Make2D(*c, "etaVsIPzMC",
2747  "#eta vs IP_{#it{z}} signal",
2748  kCyan+2,24, etaAxis, ipzAxis);
2749  fEtaVsIPzMC->SetYTitle("IP_{#it{z},gen}");
2750  fEtaVsIPzMCSel = Make2D(*c, "etaVsIPzMCSel",
2751  "#eta vs IP_{#it{z}} selected signal",
2752  kMagenta+2,24, etaAxis, ipzAxis);
2753  fEtaVsIPzMCSel->SetYTitle("IP_{#it{z},rec}");
2754  fPdgMC = Make1D(*c, "pdgMC", "Particle types",
2755  kRed+2,20, pdg);
2756  fPdgMC->GetXaxis()->LabelsOption("v");
2757  fPrimaryPdg = Make2D(*c, "primaryPdg", "Primary specie",
2758  kBlue+2, 20, pdg, deltaAxis);
2759  fSecondaryPdg = Make2D(*c, "secondaryPdg", "Secondary specie",
2760  kGreen+2, 21, pdg, deltaAxis);
2761  fPrimaryParentPdg = Make2D(*c, "primaryParentPdg",
2762  "Parent of primary specie",
2763  kBlue+2, 20, pdg, deltaAxis);
2764  fSecondaryParentPdg = Make2D(*c, "secondaryParentPdg",
2765  "Parent of secondary specie",
2766  kGreen+2, 21, pdg, deltaAxis);
2767  fPrimaryPdg->GetXaxis()->LabelsOption("v");
2768  fSecondaryPdg->GetXaxis()->LabelsOption("v");
2769  fPrimaryParentPdg->GetXaxis()->LabelsOption("v");
2770  fSecondaryParentPdg->GetXaxis()->LabelsOption("v");
2771  TAxis dIPz(ipzAxis);
2772  ScaleAxis(dIPz, .1);
2773  dIPz.SetTitle("IP_{#it{z}} - IP_{#it{z},gen}");
2774  fIPzVsMC = Make2D(*c,"ipzVsGenIPz", "Resolution vs IP_{#it{z}}",
2775  kCyan+2, 24, ipzAxis, dIPz);
2776  fIPzGen = Make1D(*c,"ipzGen", "Generated IP_{z}",
2777  kCyan+2, 25, ipzAxis);
2778  fIPzSel = Make1D(*c,"ipzSel", "Selected, generated IP_{#it{z}}",
2779  kCyan+2, 21, ipzAxis);
2780  fIPzGen->SetMarkerSize(1.2*fIPzGen->GetMarkerSize());
2781  fIPz->SetMarkerStyle(24);
2782  fIPz->SetMarkerColor(kCyan+2);
2783  fIPz->SetMarkerSize(1.4*fIPz->GetMarkerSize());
2784  fIPz->SetLineColor(kCyan+2);
2785 }
2786 //____________________________________________________________________
2787 Bool_t AliTrackletdNdetaMCTask::CentBin::FillSet(UShort_t set,
2788  Bool_t isSignal,
2789  Double_t ipZ,
2790  Double_t eta,
2791  Double_t dPhi,
2792  Double_t dPhiS,
2793  Double_t dTheta,
2794  Double_t dThetaX,
2795  Double_t delta,
2796  Double_t weight)
2797 {
2798  if (AliTrackletdNdetaTask::CentBin::FillSet(set, isSignal, ipZ, eta,
2799  dPhi, dPhiS, dTheta, dThetaX,
2800  delta, weight)) return true;
2801  HistoSet* hset = 0;
2802  if (set == kPrimaries) hset = fPrimSet;
2803  else if (set == kSecondaries) hset = fSecSet;
2804  else if (set == kCombinatorics) hset = fCombSet;
2805  else if (set == kUncorrCombinatorics) hset = fCombUSet;
2806 
2807  if (!hset) return false;
2808 
2809  hset->Fill(isSignal ? ipZ : -1000,
2810  eta, dPhi, dPhiS, dTheta, dThetaX, delta, weight);
2811 
2812  return true;
2813 }
2814 //____________________________________________________________________
2815 void AliTrackletdNdetaMCTask::CentBin::FillSpecie(Bool_t primary,
2816  Int_t pdg,
2817  Int_t parentPdg,
2818  Double_t delta,
2819  Double_t weight)
2820 {
2821  Int_t bin = FindPdg(pdg);
2822  Int_t parentBin = FindPdg(parentPdg);
2823  TH2* pdgHist = (primary ? fPrimaryPdg : fSecondaryPdg);
2824  TH2* parHist = (primary ? fPrimaryParentPdg : fSecondaryParentPdg);
2825  if (pdg >= 0) pdgHist->Fill(bin, delta, weight);
2826  if (parentPdg >= 0) parHist->Fill(parentBin, delta, weight);
2827 }
2828 //____________________________________________________________________
2829 void AliTrackletdNdetaMCTask::CentBin::FillMCIPz(Double_t ipZ, Double_t genIPz)
2830 {
2831  fIPzGen->Fill(genIPz);
2832  if (ipZ < -999) return;
2833  fIPzSel->Fill(genIPz);
2834  fIPzVsMC->Fill(ipZ, ipZ-genIPz);
2835 }
2836 //____________________________________________________________________
2837 void AliTrackletdNdetaMCTask::CentBin::FillPrimary(Double_t ipZ,
2838  Double_t genIPz,
2839  Double_t eta,
2840  Bool_t charged,
2841  Int_t pdg,
2842  Double_t weight)
2843 {
2844  // Fill PDG histogram here
2845  fPdgMC->Fill(FindPdg(pdg));
2846  if (!charged) return;
2847  fEtaVsIPzMC->Fill(eta, genIPz, weight);
2848  if (ipZ > -999) fEtaVsIPzMCSel->Fill(eta, ipZ, weight);
2849 }
2850 
2851 //____________________________________________________________________
2852 TH2*
2853 AliTrackletdNdetaMCTask::CentBin::MakeAlpha(const char* name,
2854  const char* title,
2855  Container* result,
2856  TH2* primaries) const
2857 {
2858  if (!primaries) {
2859  AliWarning("No primaries supplied for alpha calculation");
2860  return 0;
2861  }
2862  if (!result) {
2863  AliWarning("No data supplied for alpha calculations");
2864  return 0;
2865  }
2866  TH2* data = GetH2(result, "signalEst");
2867  TH2* alpha = static_cast<TH2*>(CloneAndAdd(result,primaries));
2868  alpha->SetName(name);
2869  alpha->SetTitle(title);
2870  alpha->Divide(data);
2871  return alpha;
2872 }
2873 
2874 //____________________________________________________________________
2875 TH2*
2876 AliTrackletdNdetaMCTask::CentBin::MakeAlphaMask(TH2* a1,
2877  TH2* a2,
2878  Double_t min,
2879  Double_t max) const
2880 {
2881  if (!a1 || !a2) return 0;
2882  TH2* mst = (a1 ? a1 : a2);
2883  TH2* ret = static_cast<TH2*>(mst->Clone("fiducial"));
2884  ret->SetTitle("Fiducial cut");
2885  ret->SetDirectory(0);
2886  ret->SetMarkerColor(kBlack);
2887  ret->SetMarkerSize(1);
2888  ret->SetMarkerStyle(1);
2889  ret->SetFillColor(kBlack);
2890  ret->SetFillStyle(1001);
2891  ret->SetLineColor(kBlack);
2892  ret->Reset();
2893 
2894  for (Int_t ix = 1; ix <= mst->GetNbinsX(); ix++) {
2895  for (Int_t iy = 1; iy <= mst->GetNbinsY(); iy++) {
2896  Double_t c1 = (a1 ? a1->GetBinContent(ix,iy) : min);
2897  Double_t c2 = (a2 ? a2->GetBinContent(ix,iy) : min);
2898  if ((c1 > min && c1 <= max) ||
2899  (c2 > min && c2 <= max))
2900  ret->SetBinContent(ix,iy,1);
2901  }
2902  }
2903  return ret;
2904 }
2905 
2906 //____________________________________________________________________
2907 TH1*
2908 AliTrackletdNdetaMCTask::CentBin::MakePdgProj(const char* name,
2909  const char* title,
2910  TH2* orig,
2911  Double_t cut) const
2912 {
2913  if (!orig) {
2914  AliWarningF("Nothing to make %s from", name);
2915  return 0;
2916  }
2917  Int_t first = 1;
2918  Int_t last = (cut < 0 ?
2919  orig->GetYaxis()->GetNbins() :
2920  orig->GetYaxis()->FindBin(cut-1e-6));
2921  Double_t intg = orig->Integral();
2922  TH1* proj = orig->ProjectionX(name,first,last,"e");
2923  if (!proj) {
2924  AliWarningF("Failed to project %s between 0 and %f",
2925  orig->GetName(), cut);
2926  return 0;
2927  }
2928  if (cut > 0) {
2929  proj->SetMarkerStyle(proj->GetMarkerStyle()+4);
2930  proj->SetMarkerSize(1.2*proj->GetMarkerSize());
2931  }
2932  proj->SetFillStyle(0);
2933  proj->SetFillColor(0);
2934  proj->SetTitle(title);
2935  proj->SetDirectory(0);
2936  // Scale to integral of whole thing. In this way, we can see how
2937  // many particles of each kind (primary, secondary) we cut away with
2938  // our delta cut, or in case of parents, which kind of parents are
2939  // most likely to result in ignored particles. This differs a
2940  // little from how it was done in Robertos script, where the
2941  // normalisation is done to the total within the cut. In that case,
2942  // the plot shows merely the relative abundance of particles before
2943  // and after applying the cut - not which were cut away. On the
2944  // other hand, in Robertos task the histogram filled was PDG vs
2945  // centrality - not PDG vs Delta - as is done here and what seems to
2946  // be assumed in Roberto's script.
2947  if (intg > 0) proj->Scale(1./intg);
2948  return proj;
2949 }
2950 
2951 
2952 //____________________________________________________________________
2953 AliTrackletdNdetaTask::Container*
2954 AliTrackletdNdetaMCTask::CentBin::MasterFinalize(Container* parent,
2955  TH1* ipz,
2956  Double_t deltaCut,
2957  Double_t deltaTail)
2958 {
2959  Container* result =
2960  AliTrackletdNdetaTask::CentBin::MasterFinalize(parent, ipz,
2961  deltaCut, deltaTail);
2962  TH1* ipzRec = GetH1(fContainer, "ipz"); // Original - not scaled
2963  TH1* ipzGen = GetH1(fContainer, "ipzGen");
2964 
2965  fEtaVsIPzMC = ScaleToIPz(GetH2(fContainer,"etaVsIPzMC"), ipzGen);
2966  fEtaVsIPzMCSel = ScaleToIPz(GetH2(fContainer,"etaVsIPzMCSel"),ipzRec);
2967  fIPzVsMC = GetH2(fContainer, "ipzVsGenIPz");
2968 
2969  TProfile* ipzVsMC = fIPzVsMC->ProfileX(fIPzVsMC->GetName());
2970  ipzVsMC->SetYTitle(Form("#LT%s#GT",
2971  fIPzVsMC->GetYaxis()->GetTitle()));
2972  result->Add(ipzVsMC);
2973  result->Add(fEtaVsIPzMC);
2974  result->Add(fEtaVsIPzMCSel);
2975 
2976  TH1* dNdetaAll = AverageOverIPz(fEtaVsIPzMC, "truedNdeta", 1, ipzGen);
2977  TH1* dNdetaSel = AverageOverIPz(fEtaVsIPzMCSel, "truedNdetaSel", 1, ipzRec);
2978  dNdetaAll->SetMarkerSize(1.4*dNdetaAll->GetMarkerSize());
2979  dNdetaAll->SetTitle("All events");
2980  dNdetaSel->SetTitle("Selected events");
2981  dNdetaAll->SetYTitle("d#it{N}_{ch}/d#eta");
2982  dNdetaSel->SetYTitle("d#it{N}_{ch}/d#eta");
2983  result->Add(dNdetaAll);
2984  result->Add(dNdetaSel);
2985 
2986  result->Add(MakePdgProj("allPrimaryPdg", "All primary particle PDGs",
2987  GetH2(fContainer, "primaryPdg")));
2988  result->Add(MakePdgProj("allSecondaryPdg","All secondary particle PDGs",
2989  GetH2(fContainer, "secondaryPdg")));
2990  result->Add(MakePdgProj("cutPrimaryPdg", "Selected primary particle PDGs",
2991  GetH2(fContainer, "primaryPdg"),deltaCut));
2992  result->Add(MakePdgProj("cutSecondaryPdg","Selected secondary particle PDGs",
2993  GetH2(fContainer, "secondaryPdg"),deltaCut));
2994  result->Add(MakePdgProj("allPrimaryParentPdg",
2995  "All primary particle parents PDGs",
2996  GetH2(fContainer, "primaryParentPdg")));
2997  result->Add(MakePdgProj("allSecondaryParentPdg",
2998  "All secondary particle parents PDGs",
2999  GetH2(fContainer, "secondaryParentPdg")));
3000  result->Add(MakePdgProj("cutPrimaryParentPdg",
3001  "Selected primary particle parents PDGs",
3002  GetH2(fContainer, "primaryParentPdg"),deltaCut));
3003  result->Add(MakePdgProj("cutSecondaryParentPdg",
3004  "Selected secondary particle parents PDGs",
3005  GetH2(fContainer, "secondaryParentPdg"),deltaCut));
3006 
3007 
3008 
3009  fIPzGen = static_cast<TH1*>(CloneAndAdd(result,ipzGen));
3010  fIPzSel = static_cast<TH1*>(CloneAndAdd(result,GetH1(fContainer, "ipzSel")));
3011  TGraphAsymmErrors* ipzEff = new TGraphAsymmErrors;
3012  ipzEff->SetName("ipzEff");
3013  ipzEff->SetTitle("IP_{z} efficiency");
3014  ipzEff->BayesDivide(fIPzSel, fIPzGen);
3015  ipzEff->SetMarkerStyle(20);
3016  ipzEff->SetMarkerColor(kRed+2);
3017  ipzEff->SetLineColor(kRed+2);
3018  ipzEff->SetFillColor(kRed-2);
3019  ipzEff->SetFillStyle(1001);
3020  fIPzGen->Scale(1./fIPzGen->GetEntries());
3021  fIPzSel->Scale(1./fIPzGen->GetEntries());
3022  result->Add(ipzEff);
3023 
3024  const char* subs[] = { "injection",
3025  "combinatorics",
3026  "uncorrelatedCombinatorics",
3027  0 };
3028  const char** ptr = subs;
3029  while (*ptr) {
3030  Container* sub = GetC(result, *ptr);
3031  if (!sub) {
3032  AliWarningF("Container %s not found in %s", *ptr, result->GetName());
3033  ptr++;
3034  }
3035  TH2* a1 = MakeAlpha("alpha", "#alpha all MC events",
3036  sub, fEtaVsIPzMC);
3037  TH2* a2 = MakeAlpha("alphaSel","#alpha selected MC events",sub,
3038  fEtaVsIPzMCSel);
3039  TH2* am = MakeAlphaMask(a1, a2);
3040  sub->Add(am);
3041  ptr++;
3042  }
3043  if (lDebug > 5) {
3044  result->ls();
3045  }
3046  return result;
3047 }
3048 
3049 //====================================================================
3050 void AliTrackletdNdetaMCTask::WorkerInit()
3051 {
3052  AliTrackletdNdetaTask::WorkerInit();
3053 
3054  fMCStatus = new TH1F("statusMC", "Status of task (MC part)",
3055  kMCCompleted, .5, kMCCompleted+.5);
3056  fMCStatus->SetMarkerColor(kCyan+2);
3057  fMCStatus->SetMarkerSize(2);
3058  fMCStatus->SetLineColor(kCyan+2);
3059  fMCStatus->SetFillColor(kCyan+2);
3060  fMCStatus->SetFillStyle(1001);
3061  fMCStatus->SetDirectory(0);
3062  fMCStatus->SetStats(0);
3063  fMCStatus->SetXTitle("Event have");
3064  fMCStatus->SetYTitle("# Events");
3065  fMCStatus->GetXaxis()->SetBinLabel(kAllMC, "Been seen");
3066  fMCStatus->GetXaxis()->SetBinLabel(kMCEvent, "MC Event");
3067  fMCStatus->GetXaxis()->SetBinLabel(kStack, "Kinematics stack");
3068  fMCStatus->GetXaxis()->SetBinLabel(kGenHeader, "Generator header");
3069  fMCStatus->GetXaxis()->SetBinLabel(kDataOK, "Been selected");
3070  fMCStatus->GetXaxis()->SetBinLabel(kMCTruth, "MC 'Truth'");
3071  fMCStatus->GetXaxis()->SetBinLabel(kReweighed, "Been reweighed");
3072  fMCStatus->GetXaxis()->SetBinLabel(kMCCompleted, "Completed for MC");
3073  fContainer->Add(fMCStatus);
3074 
3075  fStatus->SetFillColor(kCyan+2);
3076  fStatus->SetLineColor(kCyan+2);
3077  fStatus->SetFillStyle(0);
3078  fStatus->SetBarOffset(0.5);
3079 
3080  fIPz ->SetMarkerSize(1.2*fIPz->GetMarkerSize());
3081  fIPz ->SetMarkerStyle(24);
3082  fIPz ->SetMarkerColor(kCyan+2);
3083  fIPz ->SetLineColor(kCyan+2);
3084 
3085  fCent ->SetMarkerSize(1.2*fIPz->GetMarkerSize());
3086  fCent ->SetMarkerStyle(24);
3087  fCent ->SetMarkerColor(kCyan+2);
3088  fCent ->SetLineColor(kCyan+2);
3089 }
3090 //____________________________________________________________________
3091 Bool_t AliTrackletdNdetaMCTask::CheckEvent(Double_t& cent,
3092  const AliVVertex*& ip,
3093  AliMultiplicity*& mult,
3094  TTree*& clusters)
3095 {
3096  // Count all seen
3097  fMCStatus->Fill(kAllMC);
3098 
3099  // Get the MC event
3100  AliMCEvent* mcEvent = MCEvent();
3101  if (!mcEvent) {
3102  AliWarning("No MC event found");
3103  return false;
3104  }
3105  fMCStatus->Fill(kMCEvent);
3106 
3107  // Get the stack
3108  if (!mcEvent->Stack()) {
3109  AliWarning("No particle stack in MC event");
3110  return false;
3111  }
3112  fMCStatus->Fill(kStack);
3113 
3114  // Get the header and generator IP
3115  TArrayF genIP(3);
3116  AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
3117  if (!genHeader) {
3118  AliWarning("No generator header found in MC event");
3119  }
3120  else {
3121  genHeader->PrimaryVertex(genIP);
3122  fMCStatus->Fill(kGenHeader);
3123  }
3124 
3125  // Now do the same event selection as for "real" data.
3126  Bool_t dataOK = AliTrackletdNdetaTask::CheckEvent(cent,ip,mult,clusters);
3127  if (dataOK) fMCStatus->Fill(kDataOK);
3128 
3129  // We reweigh the MC particles
3130  if (fReweigh) {
3131  ReweighStack(cent);
3132  fMCStatus->Fill(kReweighed);
3133  }
3134 
3135  // Fill information on primaries
3136  if (FillPrimaries(cent, ip ? ip->GetZ() : -1000, genIP[2], dataOK))
3137  fMCStatus->Fill(kMCTruth);
3138 
3139  fMCStatus->Fill(kMCCompleted);
3140  return dataOK;
3141 }
3142 //____________________________________________________________________
3143 void AliTrackletdNdetaMCTask::MasterFinalize(Container* results)
3144 {
3145  AliTrackletdNdetaTask::MasterFinalize(results);
3146  // Make copies of histograms and store
3147  fMCStatus =
3148  static_cast<TH1*>(CloneAndAdd(results,GetH1(fContainer,"statusMC")));
3149  fMCStatus->Scale(1./fMCStatus->GetBinContent(1));
3150 }
3151 
3152 //____________________________________________________________________
3153 Double_t
3154 AliTrackletdNdetaMCTask::TrackletWeight(Int_t trackletNumber,
3155  const AliMultiplicity* mult) const
3156 {
3157  if (!fReweigh) return 1;
3158 
3159  Double_t weight = 1;
3160  Int_t label[] = { mult->GetLabel(trackletNumber, 0),
3161  mult->GetLabel(trackletNumber, 1) };
3162  Int_t nLayer = 2;
3163  if (label[0] == label[1]) nLayer = 1; // from the same track
3164  for (Int_t layer = 0; layer < nLayer; layer++) {
3165  if (label[layer] < 0)
3166  // If invalid label, give this tracklet 0 weight
3167  return 0;
3168  TParticle* particle = MCEvent()->Stack()->Particle(label[layer]);
3169  if (!particle)
3170  // If particle isn't found, give this tracklet 0 weight
3171  return 0;
3172  // Assume factorization of weights
3173  weight *= particle->GetWeight();
3174  }
3175 
3176  return weight;
3177 }
3178 //____________________________________________________________________
3179 void AliTrackletdNdetaMCTask::ReweighStack(Double_t cent) const
3180 {
3181  if (!fReweigh) return;
3182 
3183  AliStack* stack = MCEvent()->Stack();
3184  Int_t nTracks = stack->GetNtrack();
3185  for (Int_t trackNo = nTracks; trackNo--;) {
3186  TParticle* particle = stack->Particle(trackNo);
3187  if (!particle) continue;
3188  particle->SetWeight(LookupWeight(trackNo, particle, cent));
3189  }
3190 }
3191 //____________________________________________________________________
3192 Double_t AliTrackletdNdetaMCTask::LookupWeight(Int_t trackNo,
3193  TParticle* particle,
3194  Double_t cent) const
3195 {
3196  if (!particle)
3197  // If we do not have a particle, give this tracklet 0 weight!
3198  return 0;
3199  AliStack* stack = MCEvent()->Stack();
3200  if (!stack->IsPhysicalPrimary(trackNo)) {
3201  // If this is not a primary, recursively call this member
3202  // function until we do find the primary mother.
3203  Int_t mother = particle->GetFirstMother();
3204  return LookupWeight(mother, stack->Particle(mother), cent);
3205  }
3206  // Now, we can do a look-up on the primary particle properties to
3207  // get the associated weight.
3208  return LookupWeight(particle, cent);
3209 }
3210 
3211 //____________________________________________________________________
3212 Bool_t AliTrackletdNdetaMCTask::FillPrimaries(Double_t cent,
3213  Double_t ipZ,
3214  Double_t genIPz,
3215  Bool_t dataOK)
3216 {
3217  AliInfoF("Fill primaries for ipZ=%f, genIPz=%f centrality=%5.1f%% (%s)",
3218  ipZ, genIPz, cent, (dataOK ? "selected" : "ignored"));
3219  Int_t nAcc = 0;
3220  TIter next(fCentBins);
3221  CentBin* bin = 0;
3222  TList toRun;
3223  while ((bin = static_cast<CentBin*>(next()))) {
3224  if (!bin->Accept(cent)) continue; // Not in range for this bin
3225  toRun.Add(bin);
3226  bin->FillMCIPz(dataOK ? ipZ : -1000, genIPz);
3227  nAcc++;
3228  }
3229  // If there's no bins to fill, give up now
3230  if (nAcc <= 0) {
3231  AliWarningF("No centrality bins found for %f%%", cent);
3232  return false;
3233  }
3234 
3235  AliStack* stack = MCEvent()->Stack();
3236  Int_t nTracks = stack->GetNtrack();
3237  for (Int_t trackNo = nTracks; trackNo--; ) {
3238  if (!stack->IsPhysicalPrimary(trackNo)) {
3239  // AliWarningF("Track # %6d not a primary", trackNo);
3240  // Not a primary, go on
3241  continue;
3242  }
3243  // Get the particle
3244  TParticle* particle = stack->Particle(trackNo);
3245  if (!particle) {
3246  AliWarningF("No particle found for track # %d", trackNo);
3247  continue;
3248  }
3249  // Get theta
3250  Double_t theta = particle->Theta();
3251  // Check for beam-like particle
3252  if (theta < 1e-6 || TMath::Abs(theta-TMath::Pi()) < 1e-6) {
3253  AliWarningF("Track # %6d is beam-like (%f)", trackNo,
3254  TMath::RadToDeg()*theta);
3255  continue;
3256  }
3257  // Get pseudorapidity, transverse momentum, weight, and type
3258  Double_t eta = particle->Eta();
3259  Double_t pT = particle->Pt();
3260  Double_t weight = particle->GetWeight();
3261  Int_t pdgID = TMath::Abs(particle->GetPdgCode());
3262  TParticlePDG* pdg = particle->GetPDG();
3263  Bool_t charge = (pdg ? pdg->Charge() != 0 : false);
3264  // @todo Fill pT spectra of charged here
3265 
3266  TIter nextB(&toRun);
3267  while ((bin = static_cast<CentBin*>(nextB()))) {
3268  // AliInfoF("Filling track %6d into %s", trackNo, bin->Name());
3269  bin->FillPrimary(dataOK ? ipZ : -1000, genIPz,
3270  eta, charge, pdgID, weight);
3271  }
3272  }
3273  return true;
3274 }
3275 
3276 //____________________________________________________________________
3277 Bool_t AliTrackletdNdetaMCTask::HaveCommonParent(const Float_t* labels0,
3278  const Float_t* labels1) const
3279 {
3280  // Max search depth
3281  const Int_t kMaxParents = 50;
3282 
3283  // Arrays and counters of parents
3284  Int_t parents[2][kMaxParents];
3285  Int_t nParents[] = { 0, 0 };
3286 
3287  // A convinience (double) array
3288  const Float_t* labels[2] = { labels0, labels1 };
3289 
3290  // Total number of tracks
3291  Int_t nTracks = MCEvent()->Stack()->GetNtrack();
3292 
3293  // Loop over layers
3294  for (Int_t layer = 0; layer < 2; layer++) {
3295  // Loop over labels - never more than 3
3296  for (Int_t labelNo = 0; labelNo < 3; labelNo++) {
3297  Int_t label = Int_t(labels[layer][labelNo]);
3298  // Check validity
3299  if (label < 0 || label >= nTracks) continue;
3300 
3301  while (nParents[layer] < kMaxParents) {
3302  // Set parent particle
3303  parents[layer][nParents[layer]++] = label;
3304  // Get the particle
3305  TParticle* particle = MCEvent()->Stack()->Particle(label);
3306  // if there's no particle, break out
3307  if (!particle) break;
3308  // Look at mother of this particle
3309  label = particle->GetFirstMother();
3310  // If there's no mother, break out
3311  if (label < 1 || label > nTracks) break;
3312  }
3313  }
3314  }
3315  // Now compare the two arrays of parents to see if we can find a
3316  // match
3317  for (Int_t label0No = nParents[0]; label0No--; ) {
3318  for (Int_t label1No = nParents[1]; label1No--; ) {
3319  if (parents[0][label0No] == parents[1][label1No])
3320  // Found common parent
3321  return true;
3322  }
3323  }
3324  return false;
3325 }
3326 
3327 //____________________________________________________________________
3328 void AliTrackletdNdetaMCTask::FillBins(TList& toRun,
3329  AliITSMultRecBg* reco,
3330  AliMultiplicity* mult,
3331  Int_t trackletNumber,
3332  Bool_t isSignal,
3333  Double_t ipz,
3334  Double_t eta,
3335  Double_t dPhi,
3336  Double_t dPhiS,
3337  Double_t dTheta,
3338  Double_t dThetaX,
3339  Double_t delta,
3340  Double_t weight) const
3341 {
3342  UShort_t mode = FindMode(reco);
3343  if (mode >= 999) return;
3344 
3345  UShort_t which = 999;
3346  Bool_t uncorr = false;
3347  Int_t pdg = -1;
3348  Int_t parentPdg = -1;
3349  if (mode == kNormal) { // Perhaps need combinatorics from ESD
3350  // Get MC labels, and derive the type off the labels
3351  // - Labels are equal, check if primary (0) of secondary (1)
3352  // - Labels are not equal - combinatorial background (2)
3353  // Do this when looking at normal reconstruction data
3354  Int_t label0 = mult->GetLabel(trackletNumber, 0);
3355  Int_t label1 = mult->GetLabel(trackletNumber, 1);
3356  which = kCombinatorics; // Assume combinatorial background
3357  if (label0 == label1)
3358  which = (MCEvent()->Stack()->IsPhysicalPrimary(label0) ?
3359  kPrimaries: kSecondaries);
3360 
3361  if (reco && which == kCombinatorics) {
3362  // Get tracklet parameters from reconstructor
3363  Float_t* trl = reco->GetTracklet(trackletNumber);
3364  // Get cluster identifiers
3365  Int_t cl0 = Int_t(trl[AliITSMultReconstructor::kClID1]);
3366  Int_t cl1 = Int_t(trl[AliITSMultReconstructor::kClID2]);
3367  // Get array of labels
3368  Float_t* labels0 = (reco->GetClusterOfLayer(0,cl0)+
3369  AliITSMultRecBg::kClMC0);
3370  Float_t* labels1 = (reco->GetClusterOfLayer(1,cl1)+
3371  AliITSMultRecBg::kClMC0);
3372  uncorr = !HaveCommonParent(labels0, labels1);
3373  }
3374  if (which != kCombinatorics) {
3375  // Track to mother particle
3376  AliStack* stack = MCEvent()->Stack();
3377  Int_t nTracks = stack->GetNtrack();
3378  TParticle* particle = stack->Particle(label0);
3379  pdg = TMath::Abs(particle->GetPdgCode());
3380  parentPdg = -1;
3381  Int_t parentID = particle->GetFirstMother();
3382  while (parentID >= 0 && parentID < nTracks) {
3383  particle = stack->Particle(parentID);
3384  parentPdg = TMath::Abs(particle->GetPdgCode());
3385  parentID = particle->GetFirstMother();
3386  }
3387  }
3388  }
3389 
3390  // Re-implemented from normal case - except part with the
3391  // uncorrelated background.
3392  TIter nextB(&toRun);
3393  CentBin* bin = 0;
3394  while ((bin = static_cast<CentBin*>(nextB()))) {
3395  // Normal operation - as per real data
3396  bin->FillSet(mode, isSignal, ipz, eta,
3397  dPhi, dPhiS, dTheta, dThetaX, delta, weight);
3398  if (which >= 999) continue;
3399 
3400  // Additional stuff for MC - i.e., fill primary, secondary, or
3401  // combinatorial distributions
3402  bin->FillSet(which, isSignal, ipz, eta, dPhi, dPhiS,
3403  dTheta, dThetaX, delta, weight);
3404 
3405  if (which != kCombinatorics)
3406  bin->FillSpecie(which == kPrimaries, pdg, parentPdg, delta, weight);
3407  if (!uncorr) continue;
3408  bin->FillSet(kUncorrCombinatorics, isSignal, ipz, eta, dPhi, dPhiS,
3409  dTheta, dThetaX, delta, weight);
3410  }
3411 }
3412 
3413 //====================================================================
3414 Int_t* AliTrackletdNdetaMCTask::PdgArray(Int_t& size)
3415 {
3416  static Int_t codes[] = {
3417  211, // #pi^{+}
3418  2212, // p
3419  321, // K^{+}
3420  323, // K^{*+}
3421  11, // e^{-}
3422  13, // #mu^{-}
3423  213, // #rho^{+}
3424  411, // D^{+}
3425  413, // D^{*+}
3426  431, // D_{s}^{+}
3427  433, // D_{s}^{*+}
3428  1114, // #Delta^{-}
3429  2214, // #Delta^{+}
3430  2224, // #Delta^{++}
3431  3112, // #Sigma^{-}
3432  3222, // #Sigma^{+}
3433  3114, // #Sigma^{*-}
3434  3224, // #Sigma^{*+}
3435  4214, // #Sigma^{*+}_{c}
3436  4224, // #Sigma^{*++}_{c}
3437  3312, // #Xi^{-}
3438  3314, // #Xi^{*-}
3439  4122, // #Lambda^{+}_{c}
3440  2112, // n
3441  2114, // #Delta^{0}
3442  22, // #gamma
3443  310, // K^{0}_{S}
3444  130, // K^{0}_{L}
3445  311, // K^{0}
3446  313, // K^{*}
3447  221, // #eta
3448  111, // #pi^{0}
3449  113, // #rho^{0}
3450  333, // #varphi
3451  331, // #eta'
3452  223, // #omega
3453  3122, // #Lambda
3454  3212, // #Sigma^{0}
3455  4114, // #Sigma^{*0}_{c}
3456  3214, // #Sigma^{*0}
3457  421, // D^{0}
3458  423, // D^{*0}
3459  3322, // #Xi_{0}
3460  3324, // #Xi^{*0}
3461  4132, // #Xi^{0}_{c}
3462  4314, // #Xi^{*0}_{c}
3463  1000000000 // Nuclei
3464  };
3465  size = sizeof(codes) / sizeof(Int_t);
3466  return codes; // Others
3467 }
3468 
3469 //____________________________________________________________________
3470 Int_t* AliTrackletdNdetaMCTask::SortedPdgArray(Int_t& size)
3471 {
3472  static Int_t nCodes = 0;
3473  static Int_t* codes = 0;
3474  static Int_t* sorted = 0;
3475  if (sorted) {
3476  size = nCodes;
3477  return sorted;
3478  }
3479 
3480  codes = PdgArray(nCodes);
3481  size = nCodes;
3482  sorted = new Int_t[nCodes];
3483  Int_t* idx = new Int_t[nCodes];
3484  TMath::Sort(nCodes, codes, idx, false);
3485  for (Int_t i = 0; i < nCodes; i++) {
3486  sorted[i] = codes[idx[i]];
3487  }
3488  delete [] idx;
3489  return sorted;
3490 }
3491 
3492 //____________________________________________________________________
3493 Int_t AliTrackletdNdetaMCTask::FindPdg(Int_t pdg)
3494 {
3495  Int_t size = 0;
3496  Int_t* array = SortedPdgArray(size);
3497  Int_t idx = TMath::BinarySearch(size, array, pdg);
3498  // Printf("Lookup of PDG %7d -> %d", pdg, idx);
3499  if (idx == size-1) return idx;
3500  if (array[idx] != pdg) return size;
3501  return idx;
3502 }
3503 
3504 //____________________________________________________________________
3505 TAxis* AliTrackletdNdetaMCTask::MakePdgAxis()
3506 {
3507  static TAxis* axis = 0;
3508  if (axis) return axis;
3509 
3510  Int_t size = 0;
3511  Int_t* array = SortedPdgArray(size);
3512  // Printf("Got array of size %d", size);
3513  axis = new TAxis(size+1, -.5, size+.5);
3514  // axis->LabelOption("v");
3515  FixAxis(*axis, "");
3516  TDatabasePDG* pdgDb = TDatabasePDG::Instance();
3517  for (Int_t i = 1; i < size; i++) {
3518  TString name = "?";
3519  Int_t pdg = array[i-1];
3520  TParticlePDG* pdgP = pdgDb->GetParticle(pdg);
3521  if (pdgP) name = pdgP->GetName();
3522  // Printf("%3d/%3d: %6d - %s", i-1, size, pdg, name.Data());
3523  axis->SetBinLabel(i, name.Data());
3524  }
3525  axis->SetBinLabel(size, "Nuclei");
3526  axis->SetBinLabel(size+1, "Unknown");
3527  return axis;
3528 }
3529 
3530 #endif
3531 
3600 /*
3601  * Result calculation in case of use of MC label combinatorical
3602  * background (Delta distributions ignored).
3603  *
3604  * P C P C
3605  * R = ---------- (1-k--)M = ---- (1-k--)M
3606  * (1-C/M*)M* M* M*-C M*
3607  *
3608  * where
3609  *
3610  * - P is the primary particle distribution
3611  * - C is the distribution of tracklets with two mothers
3612  * - M* is the observed distribution in simulated data
3613  * - M is the observed distribution in real data
3614  * - k is a scaling constant (1.3?)
3615  *
3616  * In case k=1, then we have
3617  *
3618  * P
3619  * R = -- M
3620  * M*
3621  *
3622  * and we see that P/M* is an effective correction for secondaries and
3623  * other algorithmic effects.
3624  */
3697 // Local Variables:
3698 // mode: C++
3699 // End:
3700 //
3701 // EOF
3702 //
color
Int_t color[]
Definition: DrawAnaCaloTrackQA.C:21
charge
Int_t charge
Definition: MuonTrackingEfficiency.C:54
pdg
Int_t pdg
Definition: charmCutsOptimization.C:34
AliTrackletdNdetaUtils::CopyH2
static TH2 * CopyH2(Container *parent, const char *name, const char *newName=0)
Definition: AliTrackletdNdetaUtils.C:740
AliTrackletdNdetaTask::FillBins
virtual void FillBins(TList &toRun, AliITSMultRecBg *reco, AliMultiplicity *mult, Int_t trackletNumber, Bool_t isSignal, Double_t ipz, Double_t eta, Double_t dPhi, Double_t dPhiS, Double_t dTheta, Double_t dThetaX, Double_t delta, Double_t weight) const
Definition: AliTrackletdNdetaTask.C:1951
AliTrackletdNdetaMCTask::MakeCentBin
virtual AliTrackletdNdetaTask::CentBin * MakeCentBin(Float_t c1, Float_t c2, UShort_t recFlags)
Definition: AliTrackletdNdetaTask.C:2474
AliTrackletdNdetaTask::InitCentBins
virtual void InitCentBins(Container *existing)
Definition: AliTrackletdNdetaTask.C:1543
AliTrackletdNdetaTask::HistoSet::FinalizeInit
virtual void FinalizeInit(Container *parent)
Definition: AliTrackletdNdetaTask.C:1086
AliTrackletdNdetaTask::SetZEtaOverlapCut
void SetZEtaOverlapCut(Double_t x=0.05)
Definition: AliTrackletdNdetaTask.C:711
AliTrackletdNdetaTask::CentBin::FillIPz
virtual void FillIPz(Double_t ipZ)
Definition: AliTrackletdNdetaTask.C:363
Double_t
double Double_t
Definition: External.C:58
AliTrackletdNdetaTask::SetPhiAxis
void SetPhiAxis(Int_t n, Double_t min, Double_t max)
Definition: AliTrackletdNdetaTask.C:607
title
const char * title
Definition: MakeQAPdf.C:26
AliTrackletdNdetaMCTask::CentBin::FillMCIPz
virtual void FillMCIPz(Double_t ipZ, Double_t genIPz)
Definition: AliTrackletdNdetaTask.C:2829
AliTrackletdNdetaMCTask::FillBins
virtual void FillBins(TList &toRun, AliITSMultRecBg *reco, AliMultiplicity *mult, Int_t trackletNumber, Bool_t isSignal, Double_t ipz, Double_t eta, Double_t dPhi, Double_t dPhiS, Double_t dTheta, Double_t dThetaX, Double_t delta, Double_t weight) const
Definition: AliTrackletdNdetaTask.C:3328
AliTrackletdNdetaTask::CentBin::MasterFinalize
virtual Container * MasterFinalize(Container *parent, TH1 *ipz, Double_t deltaCut, Double_t deltaTail)
Definition: AliTrackletdNdetaTask.C:1318
AliTrackletdNdetaTask::SetDeltaCut
void SetDeltaCut(Double_t x=1.5)
Definition: AliTrackletdNdetaTask.C:675
AliTrackletdNdetaTask::CentBin::FillSet
virtual Bool_t FillSet(UShort_t set, Bool_t isSignal, Double_t ipZ, Double_t eta, Double_t dPhi, Double_t dPhiS, Double_t dTheta, Double_t dThetaX, Double_t delta, Double_t weight)
Definition: AliTrackletdNdetaTask.C:1348
AliTrackletdNdetaMCTask::TrackletWeight
virtual Double_t TrackletWeight(Int_t trackletNumber, const AliMultiplicity *mult) const
Definition: AliTrackletdNdetaTask.C:3154
AliTrackletdNdetaUtils::AverageOverIPz
static TH1 * AverageOverIPz(TH2 *h, const char *name, UShort_t mode, TH1 *ipz, TH2 *mask=0)
Definition: AliTrackletdNdetaUtils.C:994
AliTrackletdNdetaTask::SetReconstructionMode
void SetReconstructionMode(const TString &mode)
Definition: AliTrackletdNdetaTask.C:535
AliTrackletdNdetaMCTask::MasterFinalize
virtual void MasterFinalize(Container *results)
Definition: AliTrackletdNdetaTask.C:3143
AliTrackletdNdetaUtils::Make1D
static TH1 * Make1D(Container &c, const TString &name, const TString &title, Color_t color, Style_t style, const TAxis &xAxis)
Definition: AliTrackletdNdetaUtils.C:751
AliTrackletdNdetaTask::SubBase::MasterFinalize
virtual Container * MasterFinalize(Container *parent, TH1 *ipz, Double_t deltaCut, Double_t deltaTail)=0
AliTrackletdNdetaMCTask::GetCDBReferenceURL
virtual const char * GetCDBReferenceURL() const
Definition: AliTrackletdNdetaTask.C:2531
AliTrackletdNdetaMCTask::CentBin::FillPrimary
virtual void FillPrimary(Double_t ipZ, Double_t genIPz, Double_t eta, Bool_t charged, Int_t pdg, Double_t weight)
Definition: AliTrackletdNdetaTask.C:2837
AliTrackletdNdetaUtils::SetAxis
static void SetAxis(TAxis &axis, Int_t n, Double_t *borders)
Definition: AliTrackletdNdetaUtils.C:879
AliTrackletdNdetaTask::SetTailDelta
void SetTailDelta(Double_t x=5)
Definition: AliTrackletdNdetaTask.C:687
AliTrackletdNdetaUtils::RatioE
static Double_t RatioE(Double_t n, Double_t en, Double_t d, Double_t ed, Double_t &er)
Definition: AliTrackletdNdetaUtils.C:1120
AliTrackletdNdetaMCTask::SortedPdgArray
static Int_t * SortedPdgArray(Int_t &size)
Definition: AliTrackletdNdetaTask.C:3470
AliTrackletdNdetaMCTask::ClassDef
ClassDef(AliTrackletdNdetaMCTask, 1)
AliTrackletdNdetaTask::SubBase::Container
AliTrackletdNdetaTask::Container Container
Definition: AliTrackletdNdetaTask.C:77
c
TCanvas * c
Definition: TestFitELoss.C:172
AliTrackletdNdetaTask::SetPhiAxis
void SetPhiAxis(Int_t n, Double_t max)
Definition: AliTrackletdNdetaTask.C:617
AliTrackletdNdetaMCTask::CentBin::MakeAlpha
virtual TH2 * MakeAlpha(const char *name, const char *title, Container *result, TH2 *primaries) const
Definition: AliTrackletdNdetaTask.C:2853
AliTrackletdNdetaTask::SetCentralityAxis
void SetCentralityAxis(Int_t n, Double_t *bins)
Definition: AliTrackletdNdetaTask.C:562
AliTrackletdNdetaTask::SubBase::CreateContainer
virtual Container * CreateContainer(Container &parent)
Definition: AliTrackletdNdetaTask.C:1050
nCentBins
Int_t nCentBins
Definition: CorrectSpectraMultiMCBG.C:184
AliTrackletdNdetaMCTask::FindPdg
static Int_t FindPdg(Int_t pdg)
Definition: AliTrackletdNdetaTask.C:3493
AliTrackletdNdetaTask::SetEtaAxis
void SetEtaAxis(const TString &spec)
Definition: AliTrackletdNdetaTask.C:642
AliTrackletdNdetaMCTask::CheckEvent
virtual Bool_t CheckEvent(Double_t &cent, const AliVVertex *&ip, AliMultiplicity *&mult, TTree *&clusters)
Definition: AliTrackletdNdetaTask.C:3091
AliTrackletdNdetaTask::CheckEvent
virtual Bool_t CheckEvent(Double_t &cent, const AliVVertex *&ip, AliMultiplicity *&mult, TTree *&clusters)
Definition: AliTrackletdNdetaTask.C:2101
AliTrackletdNdetaUtils::Integrate
static Double_t Integrate(TH1 *h, Double_t min, Double_t max, Double_t &err)
Definition: AliTrackletdNdetaUtils.C:1097
AliTrackletdNdetaTask::SetIPzAxis
void SetIPzAxis(Int_t n, Double_t min, Double_t max)
Definition: AliTrackletdNdetaTask.C:582
AliTrackletdNdetaTask::SubBase::Name
virtual const char * Name() const =0
AliTrackletdNdetaTask::GetCDBReferenceURL
virtual const char * GetCDBReferenceURL() const
Definition: AliTrackletdNdetaTask.C:811
AliTrackletdNdetaTask::GetCDBReferenceRun
virtual Int_t GetCDBReferenceRun() const
Definition: AliTrackletdNdetaTask.C:805
AliTrackletdNdetaMCTask::CentBin::MakePdgProj
virtual TH1 * MakePdgProj(const char *name, const char *title, TH2 *orig, Double_t cut=-1) const
Definition: AliTrackletdNdetaTask.C:2908
AliTrackletdNdetaUtils::ScaleToIPz
static TH2 * ScaleToIPz(TH2 *h, TH1 *ipZ, Bool_t full=false)
Definition: AliTrackletdNdetaUtils.C:957
AliTrackletdNdetaMCTask::CentBin::MasterFinalize
virtual Container * MasterFinalize(Container *parent, TH1 *ipz, Double_t deltaCut, Double_t deltaTail)
Definition: AliTrackletdNdetaTask.C:2954
AliTrackletdNdetaUtils::PrintAxis
static void PrintAxis(const TAxis &axis, Int_t nSig=2, const char *alt=0)
Definition: AliTrackletdNdetaUtils.C:941
kData
TString kData
Declare data MC or deltaAOD.
Definition: AddTaskCaloTrackCorr.C:21
AliTrackletdNdetaMCTask::CentBin::FillSpecie
virtual void FillSpecie(Bool_t primary, Int_t pdg, Int_t parentPdg, Double_t delta, Double_t weight)
Definition: AliTrackletdNdetaTask.C:2815
AliTrackletdNdetaMCTask::FillPrimaries
virtual Bool_t FillPrimaries(Double_t cent, Double_t ipZ, Double_t genIPz, Bool_t dataOK)
Definition: AliTrackletdNdetaTask.C:3212
AliTrackletdNdetaUtils::GetC
static Container * GetC(Container *parent, const char *name)
Definition: AliTrackletdNdetaUtils.C:686
AliTrackletdNdetaUtils::CopyH1
static TH1 * CopyH1(Container *parent, const char *name, const char *newName=0)
Definition: AliTrackletdNdetaUtils.C:729
AliTrackletdNdetaTask::SetDThetaWindow
void SetDThetaWindow(Double_t x=0.025)
Definition: AliTrackletdNdetaTask.C:693
AliTrackletdNdetaTask::SetDPhiWindow
void SetDPhiWindow(Double_t x=0.06)
Definition: AliTrackletdNdetaTask.C:699
AliTrackletdNdetaTask::SetScaleDThetaCut
void SetScaleDThetaCut(Double_t x=-1)
Definition: AliTrackletdNdetaTask.C:669
AliTrackletdNdetaTask::ProcessTracklets
virtual Bool_t ProcessTracklets(TList &toRun, AliITSMultRecBg *reco, Double_t cent, Double_t ipZ, AliMultiplicity *mult)
Definition: AliTrackletdNdetaTask.C:1884
Int_t
int Int_t
Definition: External.C:63
AliTrackletdNdetaTask::FindIP
virtual const AliVVertex * FindIP(AliVEvent *event, Double_t maxDispersion=0.04, Double_t maxZError=0.25)
Definition: AliTrackletdNdetaTask.C:1833
AliTrackletdNdetaTask::MasterFinalize
virtual void MasterFinalize(Container *results)
Definition: AliTrackletdNdetaTask.C:1750
UInt_t
unsigned int UInt_t
Definition: External.C:33
AliTrackletdNdetaTask::SubBase::WorkerInit
virtual void WorkerInit(Container &parent, const TAxis &etaAxis, const TAxis &ipzAxis, const TAxis &deltaAxis, const TAxis &dThetaAxis, const TAxis &dPhiAxis)=0
AliTrackletdNdetaMCTask::CentBin::operator=
CentBin & operator=(const CentBin &)
Definition: AliTrackletdNdetaTask.C:2430
AliTrackletdNdetaTask::Reconstruct
virtual Bool_t Reconstruct(TList &toRun, UShort_t mode, TTree *clusters, Double_t cent, const AliVVertex *ip)
Definition: AliTrackletdNdetaTask.C:1975
Float_t
float Float_t
Definition: External.C:68
AliTrackletdNdetaTask::SetMaxDelta
void SetMaxDelta(Double_t x=25)
Definition: AliTrackletdNdetaTask.C:681
AliTrackletdNdetaUtils::FixAxis
static void FixAxis(TAxis &axis, const char *title=0)
Definition: AliTrackletdNdetaUtils.C:849
AliTrackletdNdetaTask::HistoSet::WorkerInit
virtual void WorkerInit(Container &parent, const TAxis &etaAxis, const TAxis &ipzAxis, const TAxis &deltaAxis, const TAxis &dThetaAxis, const TAxis &dPhiAxis)
Definition: AliTrackletdNdetaTask.C:1062
AliTrackletdNdetaUtils::ScaleAxis
static void ScaleAxis(TAxis &axis, Double_t fact=1)
Definition: AliTrackletdNdetaUtils.C:865
AliTrackletdNdetaTask::SetReconstructionMode
void SetReconstructionMode(UShort_t flags)
Definition: AliTrackletdNdetaTask.C:549
AliTrackletdNdetaTask::SetEtaAxis
void SetEtaAxis(Int_t n, Double_t min, Double_t max)
Definition: AliTrackletdNdetaTask.C:628
AliTrackletdNdetaMCTask::ReweighStack
void ReweighStack(Double_t cent) const
Definition: AliTrackletdNdetaTask.C:3179
AliTrackletdNdetaTask::FindCentrality
virtual Double_t FindCentrality(AliVEvent *event)
Definition: AliTrackletdNdetaTask.C:1810
AliTrackletdNdetaMCTask::HaveCommonParent
Bool_t HaveCommonParent(const Float_t *labels0, const Float_t *labels1) const
Definition: AliTrackletdNdetaTask.C:3277
AliTrackletdNdetaMCTask::AliTrackletdNdetaMCTask
AliTrackletdNdetaMCTask(const AliTrackletdNdetaMCTask &o)
Definition: AliTrackletdNdetaTask.C:2487
AliTrackletdNdetaMCTask::Print
virtual void Print(Option_t *option) const
Definition: AliTrackletdNdetaTask.C:2505
AliTrackletdNdetaMCTask::PdgArray
static Int_t * PdgArray(Int_t &size)
Definition: AliTrackletdNdetaTask.C:3414
AliTrackletdNdetaTask::ProcessEvent
virtual void ProcessEvent(Double_t cent, const AliVVertex *ip, AliMultiplicity *mult, TTree *clusters)
Definition: AliTrackletdNdetaTask.C:2168
AliTrackletdNdetaTask::SetIPzAxis
void SetIPzAxis(Int_t n, Double_t max)
Definition: AliTrackletdNdetaTask.C:592
AliTrackletdNdetaUtils::Scale
static TH1 * Scale(TH1 *h, Double_t x, Double_t xe)
Definition: AliTrackletdNdetaUtils.C:610
AliTrackletdNdetaTask::CentBin::EstimateBackground
virtual Container * EstimateBackground(Container *dataCon, HistoSet *set, Container *result, Double_t deltaCut, Double_t deltaTail)
Definition: AliTrackletdNdetaTask.C:1226
mode
Int_t mode
Definition: anaM.C:40
AliTrackletdNdetaTask::CentBin::operator=
CentBin & operator=(const CentBin &)
Definition: AliTrackletdNdetaTask.C:430
AliTrackletdNdetaTask::FillClusters
virtual void FillClusters(AliITSMultRecBg *reco)
Definition: AliTrackletdNdetaTask.C:2034
AliTrackletdNdetaMCTask::CentBin::FillSet
virtual Bool_t FillSet(UShort_t set, Bool_t isSignal, Double_t ipZ, Double_t eta, Double_t dPhi, Double_t dPhiS, Double_t dTheta, Double_t dThetaX, Double_t delta, Double_t weight)
Definition: AliTrackletdNdetaTask.C:2787
AliTrackletdNdetaTask::TrackletWeight
virtual Double_t TrackletWeight(Int_t trackletNumber, const AliMultiplicity *mult) const
Definition: AliTrackletdNdetaTask.C:877
AliTrackletdNdetaMCTask::CentBin::MakeAlphaMask
virtual TH2 * MakeAlphaMask(TH2 *a1, TH2 *a2, Double_t min=0, Double_t max=2.5) const
Definition: AliTrackletdNdetaTask.C:2876
AliTrackletdNdetaMCTask::LookupWeight
virtual Double_t LookupWeight(TParticle *particle, Double_t cent) const
Definition: AliTrackletdNdetaTask.C:2603
AliTrackletdNdetaUtils::GetD
static Double_t GetD(Container *parent, const char *name, Double_t def=-1)
Definition: AliTrackletdNdetaUtils.C:697
AliTrackletdNdetaTask::CentBin::Accept
virtual Bool_t Accept(Float_t value) const
Definition: AliTrackletdNdetaTask.C:354
AliTrackletdNdetaTask::operator=
AliTrackletdNdetaTask & operator=(const AliTrackletdNdetaTask &o)
Definition: AliTrackletdNdetaTask.C:774
AliTrackletdNdetaTask::HistoSet::Fill
virtual void Fill(Double_t ipZ, Double_t eta, Double_t dphi, Double_t dphis, Double_t dtheta, Double_t dthetax, Double_t delta, Double_t weight=1)
Definition: AliTrackletdNdetaTask.C:1145
AliTrackletdNdetaTask::SetCentralityAxis
void SetCentralityAxis(const TString &spec)
Definition: AliTrackletdNdetaTask.C:571
AliTrackletdNdetaTask::CentBin::FinalizeInit
virtual void FinalizeInit(Container *parent)
Definition: AliTrackletdNdetaTask.C:1215
nEvents
Float_t nEvents[nProd]
Definition: DrawProductionComparison.C:36
AliTrackletdNdetaTask::SetScaleDTheta
void SetScaleDTheta(Bool_t x=false)
Definition: AliTrackletdNdetaTask.C:651
AliTrackletdNdetaTask::FindMode
UShort_t FindMode(AliITSMultRecBg *reco) const
Definition: AliTrackletdNdetaTask.C:1937
TH2
Definition: External.C:220
AliTrackletdNdetaTask::Connect
virtual Bool_t Connect(const char *sumFile=0, const char *resFile=0)
Definition: AliTrackletdNdetaTask.C:1468
AliTrackletdNdetaTask::ClassDef
ClassDef(AliTrackletdNdetaTask, 1)
AliTrackletdNdetaTask::HistoSet::MasterFinalize
virtual Container * MasterFinalize(Container *parent, TH1 *ipz, Double_t deltaCut, Double_t deltaTail)
Definition: AliTrackletdNdetaTask.C:1097
AliTrackletdNdetaTask::SetCentralityMethod
void SetCentralityMethod(const TString &name)
Definition: AliTrackletdNdetaTask.C:555
AliTrackletdNdetaTask::HistoSet::Name
virtual const char * Name() const
Definition: AliTrackletdNdetaTask.C:202
AliTrackletdNdetaUtils::GetH1
static TH1 * GetH1(Container *parent, const char *name)
Definition: AliTrackletdNdetaUtils.C:675
AliTrackletdNdetaUtils.C
AliTrackletdNdetaTask::SetPhiOverlapCut
void SetPhiOverlapCut(Double_t x=0.005)
Definition: AliTrackletdNdetaTask.C:705
AliTrackletdNdetaTask::SetShiftedDPhiCut
void SetShiftedDPhiCut(Double_t x=-1)
Definition: AliTrackletdNdetaTask.C:663
AliTrackletdNdetaTask::Print
virtual void Print(Option_t *) const
Definition: AliTrackletdNdetaTask.C:1517
AliTrackletdNdetaTask::CentBin::Name
virtual const char * Name() const
Definition: AliTrackletdNdetaTask.C:1193
AliTrackletdNdetaTask::SetIPzAxis
void SetIPzAxis(const TString &spec)
Definition: AliTrackletdNdetaTask.C:596
AliTrackletdNdetaTask::SetDPhiShift
void SetDPhiShift(Double_t x=0.0045)
Definition: AliTrackletdNdetaTask.C:657
AliTrackletdNdetaTask::HistoSet::HistoSet
HistoSet(const char *name="", Color_t color=kBlack, Style_t style=20)
Definition: AliTrackletdNdetaTask.C:163
UShort_t
unsigned short UShort_t
Definition: External.C:28
AliTrackletdNdetaMCTask::LookupWeight
Double_t LookupWeight(Int_t trackNo, TParticle *particle, Double_t cent) const
Definition: AliTrackletdNdetaTask.C:3192
AliTrackletdNdetaMCTask::CentBin::WorkerInit
virtual void WorkerInit(Container &parent, const TAxis &etaAxis, const TAxis &ipzAxis, const TAxis &deltaAxis, const TAxis &dThetaAxis, const TAxis &dPhiAxis)
Definition: AliTrackletdNdetaTask.C:2734
Option_t
const char Option_t
Definition: External.C:48
AliTrackletdNdetaTask::MakeCentBin
virtual CentBin * MakeCentBin(Float_t c1, Float_t c2, UShort_t recFlags)
Definition: AliTrackletdNdetaTask.C:1002
AliTrackletdNdetaMCTask::operator=
AliTrackletdNdetaMCTask & operator=(const AliTrackletdNdetaMCTask &o)
Definition: AliTrackletdNdetaTask.C:2497
AliTrackletdNdetaTask::AliTrackletdNdetaTask
AliTrackletdNdetaTask(const AliTrackletdNdetaTask &o)
Definition: AliTrackletdNdetaTask.C:735
Bool_t
bool Bool_t
Definition: External.C:53
AliTrackletdNdetaTask::HistoSet::operator=
HistoSet & operator=(const HistoSet &)
Definition: AliTrackletdNdetaTask.C:192
AliTrackletdNdetaTask::SetPhiRotation
void SetPhiRotation(Double_t x=TMath::Pi())
Definition: AliTrackletdNdetaTask.C:717
AliTrackletdNdetaTask::SubBase::FinalizeInit
virtual void FinalizeInit(Container *parent)=0
AliTrackletdNdetaUtils::Make2D
static TH2 * Make2D(Container &c, const TString &name, const TString &title, Color_t color, Style_t style, const TAxis &xAxis, const TAxis &yAxis)
Definition: AliTrackletdNdetaUtils.C:794
AliTrackletdNdetaUtils::CloneAndAdd
static TObject * CloneAndAdd(Container *c, TObject *o)
Definition: AliTrackletdNdetaUtils.C:1081
AliTrackletdNdetaTask::SetEtaAxis
void SetEtaAxis(Int_t n, Double_t max)
Definition: AliTrackletdNdetaTask.C:638
AliTrackletdNdetaTask::CentBin::WorkerInit
virtual void WorkerInit(Container &parent, const TAxis &etaAxis, const TAxis &ipzAxis, const TAxis &deltaAxis, const TAxis &dThetaAxis, const TAxis &dPhiAxis)
Definition: AliTrackletdNdetaTask.C:1200
TH1
Definition: External.C:196
AliTrackletdNdetaUtils::GetH2
static TH2 * GetH2(Container *parent, const char *name)
Definition: AliTrackletdNdetaUtils.C:680
AliTrackletdNdetaMCTask::AliTrackletdNdetaMCTask
AliTrackletdNdetaMCTask(const char *name, Bool_t reweigh=false)
Definition: AliTrackletdNdetaTask.C:2460
AliTrackletdNdetaTask::FindClusters
virtual TTree * FindClusters(Bool_t needed)
Definition: AliTrackletdNdetaTask.C:1791