PWGCF/Correlations/Base/AliUEHistograms.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: AliUEHistograms.cxx 20164 2007-08-14 15:31:50Z morsch $ */
  17
  18 //
  19 //
  20 // encapsulates several AliUEHist objects for a full UE analysis plus additional control histograms
  21 //
  22 //
  23 // Author: Jan Fiete Grosse-Oetringhaus, Sara Vallero
  24
  25 #include "AliUEHistograms.h"
  26
  27 #include "AliCFContainer.h"
  28 #include "AliVParticle.h"
  29 #include "AliAODTrack.h"
  30
  31 #include "TList.h"
  32 #include "TCanvas.h"
  33 #include "TH2F.h"
  34 #include "TH1F.h"
  35 #include "TH3F.h"
  36 #include "TMath.h"
  37 #include "TLorentzVector.h"
  38
  39 ClassImp(AliUEHistograms)
  40
  41 const Int_t AliUEHistograms::fgkUEHists = 3;
  42
  43 AliUEHistograms::AliUEHistograms(const char* name, const char* histograms, const char* binning) :
  44   TNamed(name, name),
  45   fNumberDensitypT(0),
  46   fSumpT(0),
  47   fNumberDensityPhi(0),
  48   fCorrelationpT(0),
  49   fCorrelationEta(0),
  50   fCorrelationPhi(0),
  51   fCorrelationR(0),
  52   fCorrelationLeading2Phi(0),
  53   fCorrelationMultiplicity(0),
  54   fYields(0),
  55   fEventCount(0),
  56   fEventCountDifferential(0),
  57   fVertexContributors(0),
  58   fCentralityDistribution(0),
  59   fCentralityCorrelation(0),
  60   fITSClusterMap(0),
  61   fControlConvResoncances(0),
  62   fEfficiencyCorrectionTriggers(0),
  63   fEfficiencyCorrectionAssociated(0),
  64   fSelectCharge(0),
  65   fTriggerSelectCharge(0),
  66   fAssociatedSelectCharge(0),
  67   fTriggerRestrictEta(-1),
  68   fEtaOrdering(kFALSE),
  69   fCutConversions(kFALSE),
  70   fCutResonances(kFALSE),
  71   fRejectResonanceDaughters(-1),
  72   fOnlyOneEtaSide(0),
  73   fWeightPerEvent(kFALSE),
  74   fPtOrder(kTRUE),
  75   fTwoTrackCutMinRadius(0.8),
  76   fRunNumber(0),
  77   fMergeCount(1)
  78 {
  79   // Constructor
  80   //
  81   // the string histograms defines which histograms are created:
  82   //    1 = NumberDensitypT
  83   //    2 = SumpT
  84   //    3 = NumberDensityPhi
  85   //    4 = NumberDensityPhiCentrality (other multiplicity for Pb)
  86
  87   AliLog::SetClassDebugLevel("AliCFContainer", -1);
  88   AliLog::SetClassDebugLevel("AliCFGridSparse", -3);
  89
  90   fTwoTrackDistancePt[0] = 0;
  91   fTwoTrackDistancePt[1] = 0;
  92
  93   TString histogramsStr(histograms);
  94
  95   TString defaultBinningStr;
  96   defaultBinningStr = "eta: -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0\n"
  97     "p_t_assoc: 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 8.0\n"
  98     "p_t_leading: 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0, 10.5, 11.0, 11.5, 12.0, 12.5, 13.0, 13.5, 14.0, 14.5, 15.0, 15.5, 16.0, 16.5, 17.0, 17.5, 18.0, 18.5, 19.0, 19.5, 20.0, 20.5, 21.0, 21.5, 22.0, 22.5, 23.0, 23.5, 24.0, 24.5, 25.0, 25.5, 26.0, 26.5, 27.0, 27.5, 28.0, 28.5, 29.0, 29.5, 30.0, 30.5, 31.0, 31.5, 32.0, 32.5, 33.0, 33.5, 34.0, 34.5, 35.0, 35.5, 36.0, 36.5, 37.0, 37.5, 38.0, 38.5, 39.0, 39.5, 40.0, 40.5, 41.0, 41.5, 42.0, 42.5, 43.0, 43.5, 44.0, 44.5, 45.0, 45.5, 46.0, 46.5, 47.0, 47.5, 48.0, 48.5, 49.0, 49.5, 50.0\n"
  99     "p_t_leading_course: 0.5, 1.0, 2.0, 3.0, 4.0, 6.0, 8.0\n"
 100     "p_t_eff: 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0\n"
 101     "vertex_eff: -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10\n"
 102   ;
 103
 104   if (histogramsStr.Contains("4") || histogramsStr.Contains("5") || histogramsStr.Contains("6")) // Dphi Corr
 105   {
 106     if (histogramsStr.Contains("C"))
 107       defaultBinningStr += "multiplicity: 0, 20, 40, 60, 80, 100.1\n"; // course
 108     else
 109       defaultBinningStr += "multiplicity: 0, 1, 2, 3, 4, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100.1\n";
 110
 111     if (histogramsStr.Contains("5"))
 112       defaultBinningStr += "vertex: -7, -5, -3, -1, 1, 3, 5, 7\n";
 113     else
 114       defaultBinningStr += "vertex: -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10\n";
 115
 116     if (histogramsStr.Contains("R"))
 117       defaultBinningStr += "delta_phi: -1.570796, -1.483530, -1.396263, -1.308997, -1.221730, -1.134464, -1.047198, -0.959931, -0.872665, -0.785398, -0.698132, -0.610865, -0.523599, -0.436332, -0.349066, -0.261799, -0.174533, -0.087266, 0.0, 0.087266, 0.174533, 0.261799, 0.349066, 0.436332, 0.523599, 0.610865, 0.698132, 0.785398, 0.872665, 0.959931, 1.047198, 1.134464, 1.221730, 1.308997, 1.396263, 1.483530, 1.570796, 1.658063, 1.745329, 1.832596, 1.919862, 2.007129, 2.094395, 2.181662, 2.268928, 2.356194, 2.443461, 2.530727, 2.617994, 2.705260, 2.792527, 2.879793, 2.967060, 3.054326, 3.141593, 3.228859, 3.316126, 3.403392, 3.490659, 3.577925, 3.665191, 3.752458, 3.839724, 3.926991, 4.014257, 4.101524, 4.188790, 4.276057, 4.363323, 4.450590, 4.537856, 4.625123, 4.712389\n" // this binning starts at -pi/2 and is modulo 3
 118         "delta_eta: -2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0,2.1, 2.2, 2.3, 2.4\n"
 119       ;
 120     else // for TTR studies
 121       defaultBinningStr += "delta_phi: -1.570796, -1.483530, -1.396263, -1.308997, -1.221730, -1.134464, -1.047198, -0.959931, -0.872665, -0.785398, -0.698132, -0.610865, -0.523599, -0.436332, -0.349066, -0.261799, -0.174533, -0.087266, -0.043633, -0.021817, 0.0, 0.021817, 0.043633, 0.087266, 0.174533, 0.261799, 0.349066, 0.436332, 0.523599, 0.610865, 0.698132, 0.785398, 0.872665, 0.959931, 1.047198, 1.134464, 1.221730, 1.308997, 1.396263, 1.483530, 1.570796, 1.658063, 1.745329, 1.832596, 1.919862, 2.007129, 2.094395, 2.181662, 2.268928, 2.356194, 2.443461, 2.530727, 2.617994, 2.705260, 2.792527, 2.879793, 2.967060, 3.054326, 3.141593, 3.228859, 3.316126, 3.403392, 3.490659, 3.577925, 3.665191, 3.752458, 3.839724, 3.926991, 4.014257, 4.101524, 4.188790, 4.276057, 4.363323, 4.450590, 4.537856, 4.625123, 4.712389\n"
 122         "delta_eta: -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, -0.05, -0.025, 0, 0.025, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0\n"
 123       ;
 124   }
 125   else // UE
 126     defaultBinningStr += "multiplicity: -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 200.0\n"
 127       "delta_phi: -1.570796, -1.483530, -1.396263, -1.308997, -1.221730, -1.134464, -1.047198, -0.959931, -0.872665, -0.785398, -0.698132, -0.610865, -0.523599, -0.436332, -0.349066, -0.261799, -0.174533, -0.087266, 0.0, 0.087266, 0.174533, 0.261799, 0.349066, 0.436332, 0.523599, 0.610865, 0.698132, 0.785398, 0.872665, 0.959931, 1.047198, 1.134464, 1.221730, 1.308997, 1.396263, 1.483530, 1.570796, 1.658063, 1.745329, 1.832596, 1.919862, 2.007129, 2.094395, 2.181662, 2.268928, 2.356194, 2.443461, 2.530727, 2.617994, 2.705260, 2.792527, 2.879793, 2.967060, 3.054326, 3.141593, 3.228859, 3.316126, 3.403392, 3.490659, 3.577925, 3.665191, 3.752458, 3.839724, 3.926991, 4.014257, 4.101524, 4.188790, 4.276057, 4.363323, 4.450590, 4.537856, 4.625123, 4.712389\n"
 128       "vertex: -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10\n"
 129     ;
 130
 131   // combine customBinning with defaultBinningStr -> use customBinning where available and otherwise defaultBinningStr
 132   TString customBinning(binning);
 133   TString binningStr;
 134
 135   TObjArray* lines = defaultBinningStr.Tokenize("\n");
 136   for (Int_t i=0; i<lines->GetEntriesFast(); i++)
 137   {
 138     TString line(lines->At(i)->GetName());
 139     TString tag = line(0, line.Index(":")+1);
 140     if (!customBinning.BeginsWith(tag) && !customBinning.Contains(TString("\n") + tag))
 141       binningStr += line + "\n";
 142     else
 143       Printf("Using custom binning for %s", tag.Data());
 144   }
 145   delete lines;
 146   binningStr += customBinning;
 147
 148   if (histogramsStr.Contains("1"))
 149     fNumberDensitypT = new AliUEHist("NumberDensitypT", binningStr);
 150   if (histogramsStr.Contains("2"))
 151     fSumpT = new AliUEHist("SumpT", binningStr);
 152
 153   if (histogramsStr.Contains("3"))
 154     fNumberDensityPhi = new AliUEHist("NumberDensityPhi", binningStr);
 155   else if (histogramsStr.Contains("4"))
 156     fNumberDensityPhi = new AliUEHist("NumberDensityPhiCentrality", binningStr);
 157   else if (histogramsStr.Contains("5") || histogramsStr.Contains("6"))
 158     fNumberDensityPhi = new AliUEHist("NumberDensityPhiCentralityVtx", binningStr);
 159
 160   // do not add this hists to the directory
 161   Bool_t oldStatus = TH1::AddDirectoryStatus();
 162   TH1::AddDirectory(kFALSE);
 163
 164   if (!histogramsStr.Contains("4") && !histogramsStr.Contains("5") && !histogramsStr.Contains("6"))
 165   {
 166     fCorrelationpT  = new TH2F("fCorrelationpT", ";p_{T,lead} (MC);p_{T,lead} (RECO)", 200, 0, 50, 200, 0, 50);
 167     fCorrelationEta = new TH2F("fCorrelationEta", ";#eta_{lead} (MC);#eta_{T,lead} (RECO)", 200, -1, 1, 200, -1, 1);
 168     fCorrelationPhi = new TH2F("fCorrelationPhi", ";#varphi_{lead} (MC);#varphi_{T,lead} (RECO)", 200, 0, TMath::TwoPi(), 200, 0, TMath::TwoPi());
 169   }
 170   else
 171   {
 172     fCorrelationpT  = new TH2F("fCorrelationpT", ";Centrality;p_{T} (RECO)", 100, 0, 100.001, 200, 0, 50);
 173     fCorrelationEta = new TH2F("fCorrelationEta", ";Centrality;#eta (RECO)", 100, 0, 100.001, 200, -5, 5);
 174     fCorrelationPhi = new TH2F("fCorrelationPhi", ";Centrality;#varphi (RECO)", 100, 0, 100.001, 200, 0, TMath::TwoPi());
 175   }
 176
 177   fCorrelationR =   new TH2F("fCorrelationR", ";R;p_{T,lead} (MC)", 200, 0, 2, 200, 0, 50);
 178   fCorrelationLeading2Phi = new TH2F("fCorrelationLeading2Phi", ";#Delta #varphi;p_{T,lead} (MC)", 200, -TMath::Pi(), TMath::Pi(), 200, 0, 50);
 179   fCorrelationMultiplicity = new TH2F("fCorrelationMultiplicity", ";MC tracks;Reco tracks", 100, -0.5, 99.5, 100, -0.5, 99.5);
 180   fYields = new TH3F("fYields", ";centrality;pT;eta", 100, 0, 100, 40, 0, 20, 100, -1, 1);
 181
 182   if (!histogramsStr.Contains("4") && !histogramsStr.Contains("5") && !histogramsStr.Contains("6"))
 183   {
 184     fEventCount = new TH2F("fEventCount", ";step;event type;count", AliUEHist::fgkCFSteps+2, -2.5, -0.5 + AliUEHist::fgkCFSteps, 3, -0.5, 2.5);
 185     fEventCount->GetYaxis()->SetBinLabel(1, "ND");
 186     fEventCount->GetYaxis()->SetBinLabel(2, "SD");
 187     fEventCount->GetYaxis()->SetBinLabel(3, "DD");
 188   }
 189   else
 190   {
 191     fEventCount = new TH2F("fEventCount", ";step;centrality;count", AliUEHist::fgkCFSteps+2, -2.5, -0.5 + AliUEHist::fgkCFSteps, fNumberDensityPhi->GetEventHist()->GetNBins(1), fNumberDensityPhi->GetEventHist()->GetAxis(1, 0)->GetXbins()->GetArray());
 192   }
 193
 194   fEventCountDifferential = new TH3F("fEventCountDifferential", ";p_{T,lead};step;event type", 100, 0, 50, AliUEHist::fgkCFSteps, -0.5, -0.5 + AliUEHist::fgkCFSteps, 3, -0.5, 2.5);
 195   fEventCountDifferential->GetZaxis()->SetBinLabel(1, "ND");
 196   fEventCountDifferential->GetZaxis()->SetBinLabel(2, "SD");
 197   fEventCountDifferential->GetZaxis()->SetBinLabel(3, "DD");
 198
 199   fVertexContributors = new TH1F("fVertexContributors", ";contributors;count", 100, -0.5, 99.5);
 200
 201   if (fNumberDensityPhi)
 202   {
 203     fCentralityDistribution = new TH1F("fCentralityDistribution", ";centrality;count", fNumberDensityPhi->GetEventHist()->GetNBins(1), fNumberDensityPhi->GetEventHist()->GetAxis(1, 0)->GetXbins()->GetArray());
 204     fCentralityCorrelation = new TH2F("fCentralityCorrelation", ";centrality;multiplicity", 404, 0, 101, 200, 0, 4000);
 205   }
 206
 207   fITSClusterMap = new TH3F("fITSClusterMap", "; its cluster map; centrality; pT", 256, -0.5, 255.5, 20, 0, 100.001, 100, 0, 20);
 208
 209   fControlConvResoncances = new TH2F("fControlConvResoncances", ";id;delta mass", 3, -0.5, 2.5, 100, -0.1, 0.1);
 210
 211   TH1::AddDirectory(oldStatus);
 212 }
 213
 214 //_____________________________________________________________________________
 215 AliUEHistograms::AliUEHistograms(const AliUEHistograms &c) :
 216   TNamed(fName, fTitle),
 217   fNumberDensitypT(0),
 218   fSumpT(0),
 219   fNumberDensityPhi(0),
 220   fCorrelationpT(0),
 221   fCorrelationEta(0),
 222   fCorrelationPhi(0),
 223   fCorrelationR(0),
 224   fCorrelationLeading2Phi(0),
 225   fCorrelationMultiplicity(0),
 226   fYields(0),
 227   fEventCount(0),
 228   fEventCountDifferential(0),
 229   fVertexContributors(0),
 230   fCentralityDistribution(0),
 231   fCentralityCorrelation(0),
 232   fITSClusterMap(0),
 233   fControlConvResoncances(0),
 234   fEfficiencyCorrectionTriggers(0),
 235   fEfficiencyCorrectionAssociated(0),
 236   fSelectCharge(0),
 237   fTriggerSelectCharge(0),
 238   fAssociatedSelectCharge(0),
 239   fTriggerRestrictEta(-1),
 240   fEtaOrdering(kFALSE),
 241   fCutConversions(kFALSE),
 242   fCutResonances(kFALSE),
 243   fRejectResonanceDaughters(-1),
 244   fOnlyOneEtaSide(0),
 245   fWeightPerEvent(kFALSE),
 246   fPtOrder(kTRUE),
 247   fTwoTrackCutMinRadius(0.8),
 248   fRunNumber(0),
 249   fMergeCount(1)
 250 {
 251   //
 252   // AliUEHistograms copy constructor
 253   //
 254
 255   fTwoTrackDistancePt[0] = 0;
 256   fTwoTrackDistancePt[1] = 0;
 257
 258   ((AliUEHistograms &) c).Copy(*this);
 259 }
 260
 261 //____________________________________________________________________
 262 AliUEHistograms::~AliUEHistograms()
 263 {
 264   // Destructor
 265
 266   DeleteContainers();
 267 }
 268
 269 void AliUEHistograms::DeleteContainers()
 270 {
 271   if (fNumberDensitypT)
 272   {
 273     delete fNumberDensitypT;
 274     fNumberDensitypT = 0;
 275   }
 276
 277   if (fSumpT)
 278   {
 279     delete fSumpT;
 280     fSumpT = 0;
 281   }
 282
 283   if (fNumberDensityPhi)
 284   {
 285     delete fNumberDensityPhi;
 286     fNumberDensityPhi = 0;
 287   }
 288
 289   if (fCorrelationpT)
 290   {
 291     delete fCorrelationpT;
 292     fCorrelationpT = 0;
 293   }
 294
 295   if (fCorrelationEta)
 296   {
 297     delete fCorrelationEta;
 298     fCorrelationEta = 0;
 299   }
 300
 301   if (fCorrelationPhi)
 302   {
 303     delete fCorrelationPhi;
 304     fCorrelationPhi = 0;
 305   }
 306
 307   if (fCorrelationR)
 308   {
 309     delete fCorrelationR;
 310     fCorrelationR = 0;
 311   }
 312
 313   if (fCorrelationLeading2Phi)
 314   {
 315     delete fCorrelationLeading2Phi;
 316     fCorrelationLeading2Phi = 0;
 317   }
 318
 319   if (fCorrelationMultiplicity)
 320   {
 321     delete fCorrelationMultiplicity;
 322     fCorrelationMultiplicity = 0;
 323   }
 324
 325   if (fYields)
 326   {
 327     delete fYields;
 328     fYields = 0;
 329   }
 330
 331   if (fEventCount)
 332   {
 333     delete fEventCount;
 334     fEventCount = 0;
 335   }
 336
 337   if (fEventCountDifferential)
 338   {
 339     delete fEventCountDifferential;
 340     fEventCountDifferential = 0;
 341   }
 342
 343   if (fVertexContributors)
 344   {
 345     delete fVertexContributors;
 346     fVertexContributors = 0;
 347   }
 348
 349   if (fCentralityDistribution)
 350   {
 351     delete fCentralityDistribution;
 352     fCentralityDistribution = 0;
 353   }
 354
 355   if (fCentralityCorrelation)
 356   {
 357     delete fCentralityCorrelation;
 358     fCentralityCorrelation = 0;
 359   }
 360
 361   if (fITSClusterMap)
 362   {
 363     delete fITSClusterMap;
 364     fITSClusterMap = 0;
 365   }
 366
 367   for (Int_t i=0; i<2; i++)
 368     if (fTwoTrackDistancePt[i])
 369     {
 370       delete fTwoTrackDistancePt[i];
 371       fTwoTrackDistancePt[i] = 0;
 372     }
 373
 374   if (fControlConvResoncances)
 375   {
 376     delete fControlConvResoncances;
 377     fControlConvResoncances = 0;
 378   }
 379
 380   if (fEfficiencyCorrectionTriggers)
 381   {
 382     delete fEfficiencyCorrectionTriggers;
 383     fEfficiencyCorrectionTriggers = 0;
 384   }
 385
 386   if (fEfficiencyCorrectionAssociated)
 387   {
 388     delete fEfficiencyCorrectionAssociated;
 389     fEfficiencyCorrectionAssociated = 0;
 390   }
 391 }
 392
 393 AliUEHist* AliUEHistograms::GetUEHist(Int_t id)
 394 {
 395   // returns AliUEHist object, useful for loops
 396
 397   switch (id)
 398   {
 399     case 0: return fNumberDensitypT; break;
 400     case 1: return fSumpT; break;
 401     case 2: return fNumberDensityPhi; break;
 402   }
 403
 404   return 0;
 405 }
 406
 407 //____________________________________________________________________
 408 Int_t AliUEHistograms::CountParticles(TList* list, Float_t ptMin)
 409 {
 410   // counts the number of particles in the list with a pT above ptMin
 411   // TODO eta cut needed here?
 412
 413   Int_t count = 0;
 414   for (Int_t j=0; j<list->GetEntries(); j++)
 415     if (((AliVParticle*) list->At(j))->Pt() > ptMin)
 416       count++;
 417
 418   return count;
 419 }
 420
 421 //____________________________________________________________________
 422 void AliUEHistograms::Fill(Int_t eventType, Float_t zVtx, AliUEHist::CFStep step, AliVParticle* leading, TList* toward, TList* away, TList* min, TList* max)
 423 {
 424   // fills the UE event histograms
 425   //
 426   // this function needs the leading (track or jet or ...) and four lists of AliVParticles which contain the particles/tracks to be filled in the four regions
 427
 428   // if leading is not set, just fill event statistics
 429   if (leading)
 430   {
 431     Int_t multiplicity = 0;
 432
 433     // TODO configurable?
 434     Float_t ptMin = 0.15;
 435     if (leading->Pt() > ptMin)
 436       multiplicity++;
 437
 438     multiplicity += CountParticles(toward, ptMin);
 439     multiplicity += CountParticles(away, ptMin);
 440     multiplicity += CountParticles(min, ptMin);
 441     multiplicity += CountParticles(max, ptMin);
 442
 443     FillRegion(AliUEHist::kToward, zVtx, step, leading, toward, multiplicity);
 444     FillRegion(AliUEHist::kAway,   zVtx, step, leading, away, multiplicity);
 445     FillRegion(AliUEHist::kMin,    zVtx, step, leading, min, multiplicity);
 446     FillRegion(AliUEHist::kMax,    zVtx, step, leading, max, multiplicity);
 447
 448     Double_t vars[3];
 449     vars[0] = leading->Pt();
 450     vars[1] = multiplicity;
 451     vars[2] = zVtx;
 452     for (Int_t i=0; i<fgkUEHists; i++)
 453       if (GetUEHist(i))
 454         GetUEHist(i)->GetEventHist()->Fill(vars, step);
 455
 456     fEventCountDifferential->Fill(leading->Pt(), step, eventType);
 457   }
 458
 459   FillEvent(eventType, step);
 460 }
 461
 462 //____________________________________________________________________
 463 void AliUEHistograms::FillRegion(AliUEHist::Region region, Float_t zVtx, AliUEHist::CFStep step, AliVParticle* leading, TList* list, Int_t multiplicity)
 464 {
 465   // loops over AliVParticles in list and fills the given region at the given step
 466   //
 467   // See also Fill(...)
 468
 469   for (Int_t i=0; i<list->GetEntries(); i++)
 470   {
 471     AliVParticle* particle = (AliVParticle*) list->At(i);
 472
 473     Double_t vars[6];
 474     vars[0] = particle->Eta();
 475     vars[1] = particle->Pt();
 476     vars[2] = leading->Pt();
 477     vars[3] = multiplicity;
 478     vars[4] = leading->Phi() - particle->Phi();
 479     if (vars[4] > 1.5 * TMath::Pi())
 480       vars[4] -= TMath::TwoPi();
 481     if (vars[4] < -0.5 * TMath::Pi())
 482       vars[4] += TMath::TwoPi();
 483     vars[5] = zVtx;
 484
 485     if (fNumberDensitypT)
 486       fNumberDensitypT->GetTrackHist(region)->Fill(vars, step);
 487
 488     if (fSumpT)
 489       fSumpT->GetTrackHist(region)->Fill(vars, step, particle->Pt());
 490
 491     // fill all in toward region (is anyway as function of delta phi!)
 492     if (fNumberDensityPhi)
 493       fNumberDensityPhi->GetTrackHist(AliUEHist::kToward)->Fill(vars, step);
 494   }
 495 }
 496
 497 //____________________________________________________________________
 498 void AliUEHistograms::Fill(AliVParticle* leadingMC, AliVParticle* leadingReco)
 499 {
 500   // fills the correlation histograms
 501
 502   if (leadingMC)
 503   {
 504     fCorrelationpT->Fill(leadingMC->Pt(), leadingReco->Pt());
 505     if (leadingMC->Pt() > 0.5)
 506     {
 507       fCorrelationEta->Fill(leadingMC->Eta(), leadingReco->Eta());
 508       fCorrelationPhi->Fill(leadingMC->Phi(), leadingReco->Phi());
 509     }
 510
 511     Float_t phiDiff = leadingMC->Phi() - leadingReco->Phi();
 512     if (phiDiff > TMath::Pi())
 513       phiDiff -= TMath::TwoPi();
 514     if (phiDiff < -TMath::Pi())
 515       phiDiff += TMath::TwoPi();
 516
 517     Float_t etaDiff = leadingMC->Eta() - leadingReco->Eta();
 518
 519     fCorrelationR->Fill(TMath::Sqrt(phiDiff * phiDiff + etaDiff * etaDiff), leadingMC->Pt());
 520     fCorrelationLeading2Phi->Fill(phiDiff, leadingMC->Pt());
 521   }
 522   else
 523   {
 524     fCorrelationpT->Fill(1.0, leadingReco->Pt());
 525     if (leadingReco->Pt() > 0.5)
 526     {
 527       fCorrelationEta->Fill(0.0, leadingReco->Eta());
 528       fCorrelationPhi->Fill(0.0, leadingReco->Phi());
 529     }
 530   }
 531 }
 532
 533 //____________________________________________________________________
 534 void AliUEHistograms::FillCorrelations(Double_t centrality, Float_t zVtx, AliUEHist::CFStep step, TObjArray* particles, TObjArray* mixed, Float_t weight, Bool_t firstTime, Bool_t twoTrackEfficiencyCut, Float_t bSign, Float_t twoTrackEfficiencyCutValue, Bool_t applyEfficiency)
 535 {
 536   // fills the fNumberDensityPhi histogram
 537   //
 538   // this function need a list of AliVParticles which contain the particles/tracks to be filled
 539   //
 540   // if mixed is non-0, mixed events are filled, the trigger particle is from particles, the associated from mixed
 541   // if weight < 0, then the pt of the associated particle is filled as weight
 542
 543   Bool_t fillpT = kFALSE;
 544   if (weight < 0)
 545     fillpT = kTRUE;
 546
 547   if (twoTrackEfficiencyCut && !fTwoTrackDistancePt[0])
 548   {
 549     // do not add this hists to the directory
 550     Bool_t oldStatus = TH1::AddDirectoryStatus();
 551     TH1::AddDirectory(kFALSE);
 552
 553     fTwoTrackDistancePt[0] = new TH3F("fTwoTrackDistancePt[0]", ";#Delta#eta;#Delta#varphi^{*}_{min};#Delta p_{T}", 100, -0.15, 0.15, 100, -0.05, 0.05, 20, 0, 10);
 554     fTwoTrackDistancePt[1] = (TH3F*) fTwoTrackDistancePt[0]->Clone("fTwoTrackDistancePt[1]");
 555
 556     TH1::AddDirectory(oldStatus);
 557   }
 558
 559   // Eta() is extremely time consuming, therefore cache it for the inner loop here:
 560   TObjArray* input = (mixed) ? mixed : particles;
 561   TArrayF eta(input->GetEntriesFast());
 562   for (Int_t i=0; i<input->GetEntriesFast(); i++)
 563     eta[i] = ((AliVParticle*) input->UncheckedAt(i))->Eta();
 564
 565   // if particles is not set, just fill event statistics
 566   if (particles)
 567   {
 568     Int_t jMax = particles->GetEntriesFast();
 569     if (mixed)
 570       jMax = mixed->GetEntriesFast();
 571
 572     TH1* triggerWeighting = 0;
 573     if (fWeightPerEvent)
 574     {
 575       TAxis* axis = fNumberDensityPhi->GetTrackHist(AliUEHist::kToward)->GetGrid(0)->GetGrid()->GetAxis(2);
 576       triggerWeighting = new TH1F("triggerWeighting", "", axis->GetNbins(), axis->GetXbins()->GetArray());
 577
 578       for (Int_t i=0; i<particles->GetEntriesFast(); i++)
 579       {
 580         AliVParticle* triggerParticle = (AliVParticle*) particles->UncheckedAt(i);
 581
 582         // some optimization
 583         Float_t triggerEta = triggerParticle->Eta();
 584
 585         if (fTriggerRestrictEta > 0 && TMath::Abs(triggerEta) > fTriggerRestrictEta)
 586           continue;
 587
 588         if (fOnlyOneEtaSide != 0)
 589         {
 590           if (fOnlyOneEtaSide * triggerEta < 0)
 591             continue;
 592         }
 593
 594         if (fTriggerSelectCharge != 0)
 595           if (triggerParticle->Charge() * fTriggerSelectCharge < 0)
 596             continue;
 597
 598         triggerWeighting->Fill(triggerParticle->Pt());
 599       }
 600     }
 601
 602     // identify K, Lambda candidates and flag those particles
 603     // a TObject bit is used for this
 604     const UInt_t kResonanceDaughterFlag = 1 << 14;
 605     if (fRejectResonanceDaughters > 0)
 606     {
 607       Double_t resonanceMass = -1;
 608       Double_t massDaughter1 = -1;
 609       Double_t massDaughter2 = -1;
 610       const Double_t interval = 0.02;
 611
 612       switch (fRejectResonanceDaughters)
 613       {
 614         case 1: resonanceMass = 1.2; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // method test
 615         case 2: resonanceMass = 0.4976; massDaughter1 = 0.1396; massDaughter2 = massDaughter1; break; // k0
 616         case 3: resonanceMass = 1.115; massDaughter1 = 0.1396; massDaughter2 = 0.9383; break; // lambda
 617         default: AliFatal(Form("Invalid setting %d", fRejectResonanceDaughters));
 618       }
 619
 620       for (Int_t i=0; i<particles->GetEntriesFast(); i++)
 621         particles->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
 622       if (mixed)
 623         for (Int_t i=0; i<jMax; i++)
 624           mixed->UncheckedAt(i)->ResetBit(kResonanceDaughterFlag);
 625
 626       for (Int_t i=0; i<particles->GetEntriesFast(); i++)
 627       {
 628         AliVParticle* triggerParticle = (AliVParticle*) particles->UncheckedAt(i);
 629
 630         for (Int_t j=0; j<jMax; j++)
 631         {
 632           if (!mixed && i == j)
 633             continue;
 634
 635           AliVParticle* particle = 0;
 636           if (!mixed)
 637             particle = (AliVParticle*) particles->UncheckedAt(j);
 638           else
 639             particle = (AliVParticle*) mixed->UncheckedAt(j);
 640
 641           // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
 642           if (mixed && triggerParticle->IsEqual(particle))
 643             continue;
 644
 645           if (triggerParticle->Charge() * particle->Charge() > 0)
 646             continue;
 647
 648           Float_t mass = GetInvMassSquaredCheap(triggerParticle->Pt(), triggerParticle->Eta(), triggerParticle->Phi(), particle->Pt(), particle->Eta(), particle->Phi(), massDaughter1, massDaughter2);
 649
 650           if (TMath::Abs(mass - resonanceMass*resonanceMass) < interval*5)
 651           {
 652             mass = GetInvMassSquared(triggerParticle->Pt(), triggerParticle->Eta(), triggerParticle->Phi(), particle->Pt(), particle->Eta(), particle->Phi(), massDaughter1, massDaughter2);
 653
 654             if (mass > (resonanceMass-interval)*(resonanceMass-interval) && mass < (resonanceMass+interval)*(resonanceMass+interval))
 655             {
 656               triggerParticle->SetBit(kResonanceDaughterFlag);
 657               particle->SetBit(kResonanceDaughterFlag);
 658
 659 //            Printf("Flagged %d %d %f", i, j, TMath::Sqrt(mass));
 660             }
 661           }
 662         }
 663       }
 664     }
 665
 666     for (Int_t i=0; i<particles->GetEntriesFast(); i++)
 667     {
 668       AliVParticle* triggerParticle = (AliVParticle*) particles->UncheckedAt(i);
 669
 670       // some optimization
 671       Float_t triggerEta = triggerParticle->Eta();
 672
 673       if (fTriggerRestrictEta > 0 && TMath::Abs(triggerEta) > fTriggerRestrictEta)
 674         continue;
 675
 676       if (fOnlyOneEtaSide != 0)
 677       {
 678         if (fOnlyOneEtaSide * triggerEta < 0)
 679           continue;
 680       }
 681
 682       if (fTriggerSelectCharge != 0)
 683         if (triggerParticle->Charge() * fTriggerSelectCharge < 0)
 684           continue;
 685
 686       if (fRejectResonanceDaughters > 0)
 687         if (triggerParticle->TestBit(kResonanceDaughterFlag))
 688         {
 689 //        Printf("Skipped i=%d", i);
 690           continue;
 691         }
 692
 693       for (Int_t j=0; j<jMax; j++)
 694       {
 695         if (!mixed && i == j)
 696           continue;
 697
 698         AliVParticle* particle = 0;
 699         if (!mixed)
 700           particle = (AliVParticle*) particles->UncheckedAt(j);
 701         else
 702           particle = (AliVParticle*) mixed->UncheckedAt(j);
 703
 704         // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
 705         if (mixed && triggerParticle->IsEqual(particle))
 706           continue;
 707
 708         if (fPtOrder)
 709           if (particle->Pt() >= triggerParticle->Pt())
 710             continue;
 711
 712         if (fAssociatedSelectCharge != 0)
 713           if (particle->Charge() * fAssociatedSelectCharge < 0)
 714             continue;
 715
 716         if (fSelectCharge > 0)
 717         {
 718           // skip like sign
 719           if (fSelectCharge == 1 && particle->Charge() * triggerParticle->Charge() > 0)
 720             continue;
 721
 722           // skip unlike sign
 723           if (fSelectCharge == 2 && particle->Charge() * triggerParticle->Charge() < 0)
 724             continue;
 725         }
 726
 727         if (fEtaOrdering)
 728         {
 729           if (triggerEta < 0 && eta[j] < triggerEta)
 730             continue;
 731           if (triggerEta > 0 && eta[j] > triggerEta)
 732             continue;
 733         }
 734
 735         if (fRejectResonanceDaughters > 0)
 736           if (particle->TestBit(kResonanceDaughterFlag))
 737           {
 738 //          Printf("Skipped j=%d", j);
 739             continue;
 740           }
 741
 742         // conversions
 743         if (fCutConversions && particle->Charge() * triggerParticle->Charge() < 0)
 744         {
 745           Float_t mass = GetInvMassSquaredCheap(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), particle->Pt(), eta[j], particle->Phi(), 0.510e-3, 0.510e-3);
 746
 747           if (mass < 0.1)
 748           {
 749             mass = GetInvMassSquared(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), particle->Pt(), eta[j], particle->Phi(), 0.510e-3, 0.510e-3);
 750
 751             fControlConvResoncances->Fill(0.0, mass);
 752
 753             if (mass < 0.04*0.04)
 754               continue;
 755           }
 756         }
 757
 758         // K0s
 759         if (fCutResonances && particle->Charge() * triggerParticle->Charge() < 0)
 760         {
 761           Float_t mass = GetInvMassSquaredCheap(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), particle->Pt(), eta[j], particle->Phi(), 0.1396, 0.1396);
 762
 763           const Float_t kK0smass = 0.4976;
 764
 765           if (TMath::Abs(mass - kK0smass*kK0smass) < 0.1)
 766           {
 767             mass = GetInvMassSquared(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), particle->Pt(), eta[j], particle->Phi(), 0.1396, 0.1396);
 768
 769             fControlConvResoncances->Fill(1, mass - kK0smass*kK0smass);
 770
 771             if (mass > (kK0smass-0.02)*(kK0smass-0.02) && mass < (kK0smass+0.02)*(kK0smass+0.02))
 772               continue;
 773           }
 774         }
 775
 776         // Lambda
 777         if (fCutResonances && particle->Charge() * triggerParticle->Charge() < 0)
 778         {
 779           Float_t mass1 = GetInvMassSquaredCheap(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), particle->Pt(), eta[j], particle->Phi(), 0.1396, 0.9383);
 780           Float_t mass2 = GetInvMassSquaredCheap(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), particle->Pt(), eta[j], particle->Phi(), 0.9383, 0.1396);
 781
 782           const Float_t kLambdaMass = 1.115;
 783
 784           if (TMath::Abs(mass1 - kLambdaMass*kLambdaMass) < 0.1)
 785           {
 786             mass1 = GetInvMassSquared(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), particle->Pt(), eta[j], particle->Phi(), 0.1396, 0.9383);
 787
 788             fControlConvResoncances->Fill(2, mass1 - kLambdaMass*kLambdaMass);
 789
 790             if (mass1 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass1 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
 791               continue;
 792           }
 793           if (TMath::Abs(mass2 - kLambdaMass*kLambdaMass) < 0.1)
 794           {
 795             mass2 = GetInvMassSquared(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), particle->Pt(), eta[j], particle->Phi(), 0.9383, 0.1396);
 796
 797             fControlConvResoncances->Fill(2, mass2 - kLambdaMass*kLambdaMass);
 798
 799             if (mass2 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass2 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
 800               continue;
 801           }
 802         }
 803
 804         if (twoTrackEfficiencyCut)
 805         {
 806           // the variables & cuthave been developed by the HBT group
 807           // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
 808
 809           Float_t phi1 = triggerParticle->Phi();
 810           Float_t pt1 = triggerParticle->Pt();
 811           Float_t charge1 = triggerParticle->Charge();
 812
 813           Float_t phi2 = particle->Phi();
 814           Float_t pt2 = particle->Pt();
 815           Float_t charge2 = particle->Charge();
 816
 817           Float_t deta = triggerEta - eta[j];
 818
 819           // optimization
 820           if (TMath::Abs(deta) < twoTrackEfficiencyCutValue * 2.5 * 3)
 821           {
 822             // check first boundaries to see if is worth to loop and find the minimum
 823             Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, fTwoTrackCutMinRadius, bSign);
 824             Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 2.5, bSign);
 825
 826             const Float_t kLimit = twoTrackEfficiencyCutValue * 3;
 827
 828             Float_t dphistarminabs = 1e5;
 829             Float_t dphistarmin = 1e5;
 830             if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
 831             {
 832               for (Double_t rad=fTwoTrackCutMinRadius; rad<2.51; rad+=0.01)
 833               {
 834                 Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
 835
 836                 Float_t dphistarabs = TMath::Abs(dphistar);
 837
 838                 if (dphistarabs < dphistarminabs)
 839                 {
 840                   dphistarmin = dphistar;
 841                   dphistarminabs = dphistarabs;
 842                 }
 843               }
 844
 845               fTwoTrackDistancePt[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));
 846
 847               if (dphistarminabs < twoTrackEfficiencyCutValue && TMath::Abs(deta) < twoTrackEfficiencyCutValue)
 848               {
 849 //              Printf("Removed track pair %d %d with %f %f %f %f %f %f %f %f %f", i, j, deta, dphistarminabs, phi1, pt1, charge1, phi2, pt2, charge2, bSign);
 850                 continue;
 851               }
 852
 853               fTwoTrackDistancePt[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));
 854             }
 855           }
 856         }
 857
 858         Double_t vars[6];
 859         vars[0] = triggerEta - eta[j];
 860         vars[1] = particle->Pt();
 861         vars[2] = triggerParticle->Pt();
 862         vars[3] = centrality;
 863         vars[4] = triggerParticle->Phi() - particle->Phi();
 864         if (vars[4] > 1.5 * TMath::Pi())
 865           vars[4] -= TMath::TwoPi();
 866         if (vars[4] < -0.5 * TMath::Pi())
 867           vars[4] += TMath::TwoPi();
 868         vars[5] = zVtx;
 869
 870         if (fillpT)
 871           weight = particle->Pt();
 872
 873         Double_t useWeight = weight;
 874         if (applyEfficiency)
 875         {
 876           if (fEfficiencyCorrectionAssociated)
 877           {
 878             Int_t effVars[4];
 879             // associated particle
 880             effVars[0] = fEfficiencyCorrectionAssociated->GetAxis(0)->FindBin(eta[j]);
 881             effVars[1] = fEfficiencyCorrectionAssociated->GetAxis(1)->FindBin(vars[1]); //pt
 882             effVars[2] = fEfficiencyCorrectionAssociated->GetAxis(2)->FindBin(vars[3]); //centrality
 883             effVars[3] = fEfficiencyCorrectionAssociated->GetAxis(3)->FindBin(vars[5]); //zVtx
 884
 885             //    Printf("%d %d %d %d %f", effVars[0], effVars[1], effVars[2], effVars[3], fEfficiencyCorrectionAssociated->GetBinContent(effVars));
 886
 887             useWeight *= fEfficiencyCorrectionAssociated->GetBinContent(effVars);
 888           }
 889           if (fEfficiencyCorrectionTriggers)
 890           {
 891             Int_t effVars[4];
 892
 893             effVars[0] = fEfficiencyCorrectionTriggers->GetAxis(0)->FindBin(triggerEta);
 894             effVars[1] = fEfficiencyCorrectionTriggers->GetAxis(1)->FindBin(vars[2]); //pt
 895             effVars[2] = fEfficiencyCorrectionTriggers->GetAxis(2)->FindBin(vars[3]); //centrality
 896             effVars[3] = fEfficiencyCorrectionTriggers->GetAxis(3)->FindBin(vars[5]); //zVtx
 897             useWeight *= fEfficiencyCorrectionTriggers->GetBinContent(effVars);
 898           }
 899         }
 900
 901         if (fWeightPerEvent)
 902         {
 903           Int_t weightBin = triggerWeighting->GetXaxis()->FindBin(vars[2]);
 904 //        Printf("Using weight %f", triggerWeighting->GetBinContent(weightBin));
 905           useWeight /= triggerWeighting->GetBinContent(weightBin);
 906         }
 907
 908         // fill all in toward region and do not use the other regions
 909         fNumberDensityPhi->GetTrackHist(AliUEHist::kToward)->Fill(vars, step, useWeight);
 910
 911 //      Printf("%.2f %.2f --> %.2f", triggerEta, eta[j], vars[0]);
 912       }
 913
 914       if (firstTime)
 915       {
 916         // once per trigger particle
 917         Double_t vars[3];
 918         vars[0] = triggerParticle->Pt();
 919         vars[1] = centrality;
 920         vars[2] = zVtx;
 921
 922         Double_t useWeight = 1;
 923         if (fEfficiencyCorrectionTriggers && applyEfficiency)
 924         {
 925           Int_t effVars[4];
 926
 927           // trigger particle
 928           effVars[0] = fEfficiencyCorrectionTriggers->GetAxis(0)->FindBin(triggerEta);
 929           effVars[1] = fEfficiencyCorrectionTriggers->GetAxis(1)->FindBin(vars[0]); //pt
 930           effVars[2] = fEfficiencyCorrectionTriggers->GetAxis(2)->FindBin(vars[1]); //centrality
 931           effVars[3] = fEfficiencyCorrectionTriggers->GetAxis(3)->FindBin(vars[2]); //zVtx
 932           useWeight *= fEfficiencyCorrectionTriggers->GetBinContent(effVars);
 933         }
 934         if (fWeightPerEvent)
 935         {
 936           // leads effectively to a filling of one entry per filled trigger particle pT bin
 937           Int_t weightBin = triggerWeighting->GetXaxis()->FindBin(vars[0]);
 938 //        Printf("Using weight %f", triggerWeighting->GetBinContent(weightBin));
 939           useWeight /= triggerWeighting->GetBinContent(weightBin);
 940         }
 941
 942         fNumberDensityPhi->GetEventHist()->Fill(vars, step, useWeight);
 943
 944         // QA
 945         fCorrelationpT->Fill(centrality, triggerParticle->Pt());
 946         fCorrelationEta->Fill(centrality, triggerEta);
 947         fCorrelationPhi->Fill(centrality, triggerParticle->Phi());
 948         fYields->Fill(centrality, triggerParticle->Pt(), triggerEta);
 949
 950 /*        if (dynamic_cast<AliAODTrack*>(triggerParticle))
 951           fITSClusterMap->Fill(((AliAODTrack*) triggerParticle)->GetITSClusterMap(), centrality, triggerParticle->Pt());*/
 952       }
 953     }
 954
 955     if (triggerWeighting)
 956     {
 957       delete triggerWeighting;
 958       triggerWeighting = 0;
 959     }
 960   }
 961
 962   fCentralityDistribution->Fill(centrality);
 963   fCentralityCorrelation->Fill(centrality, particles->GetEntriesFast());
 964   FillEvent(centrality, step);
 965 }
 966
 967 //____________________________________________________________________
 968 void AliUEHistograms::FillTrackingEfficiency(TObjArray* mc, TObjArray* recoPrim, TObjArray* recoAll, TObjArray* recoPrimPID, TObjArray* recoAllPID, TObjArray* fake, Int_t particleType, Double_t centrality, Double_t zVtx)
 969 {
 970   // fills the tracking efficiency objects
 971   //
 972   // mc: all primary MC particles
 973   // recoPrim: reconstructed primaries (again MC particles)
 974   // recoAll: reconstructed (again MC particles)
 975   // recoPrim: reconstructed primaries with checks on PID (again MC particles)
 976   // recoAll: reconstructed with checks on PID (again MC particles)
 977   // particleType is: 0 for pion, 1 for kaon, 2 for proton, 3 for others
 978
 979   for (Int_t step=0; step<6; step++)
 980   {
 981     TObjArray* list = mc;
 982     if (step == 1)
 983       list = recoPrim;
 984     else if (step == 2)
 985       list = recoAll;
 986     else if (step == 3)
 987       list = recoPrimPID;
 988     else if (step == 4)
 989       list = recoAllPID;
 990     else if (step == 5)
 991       list = fake;
 992
 993     if (!list)
 994       continue;
 995
 996     for (Int_t i=0; i<list->GetEntriesFast(); i++)
 997     {
 998       AliVParticle* particle = (AliVParticle*) list->UncheckedAt(i);
 999       Double_t vars[5];
1000       vars[0] = particle->Eta();
1001       vars[1] = particle->Pt();
1002       vars[2] = particleType;
1003       vars[3] = centrality;
1004       vars[4] = zVtx;
1005
1006       for (Int_t j=0; j<fgkUEHists; j++)
1007         if (GetUEHist(j))
1008           GetUEHist(j)->GetTrackHistEfficiency()->Fill(vars, step);
1009     }
1010   }
1011 }
1012
1013 //____________________________________________________________________
1014 void AliUEHistograms::FillFakePt(TObjArray* fake, Double_t centrality)
1015 {
1016   TObjArray* tracksReco = (TObjArray*) fake->At(0);
1017   TObjArray* tracksMC = (TObjArray*) fake->At(1);
1018
1019   if (tracksReco->GetEntriesFast() != tracksMC->GetEntriesFast())
1020     AliFatal(Form("Inconsistent arrays: %d vs %d", tracksReco->GetEntriesFast(), tracksMC->GetEntriesFast()));
1021
1022   for (Int_t i=0; i<tracksReco->GetEntriesFast(); i++)
1023   {
1024     AliVParticle* particle1 = (AliVParticle*) tracksReco->At(i);
1025     AliVParticle* particle2 = (AliVParticle*) tracksMC->At(i);
1026     Double_t vars[3];
1027     vars[0] = particle1->Pt();
1028     vars[1] = particle2->Pt();
1029     vars[2] = centrality;
1030     for (Int_t j=0; j<fgkUEHists; j++)
1031       if (GetUEHist(j))
1032         GetUEHist(j)->GetMCRecoPtCorrelation()->Fill(vars[0],vars[1],vars[2]);
1033   }
1034 }
1035
1036 //____________________________________________________________________
1037 void AliUEHistograms::FillEvent(Int_t eventType, Int_t step)
1038 {
1039   // fills the number of events at the given step and the given enty type
1040   //
1041   // WARNING: This function is called from Fill, so only call it for steps where Fill is not called
1042
1043   fEventCount->Fill(step, eventType);
1044 }
1045
1046 //____________________________________________________________________
1047 void AliUEHistograms::FillEvent(Double_t centrality, Int_t step)
1048 {
1049   // fills the number of events at the given step and the given centrality
1050   //
1051   // WARNING: This function is called from Fill, so only call it for steps where Fill is not called
1052
1053   fEventCount->Fill(step, centrality);
1054 }
1055
1056 //____________________________________________________________________
1057 void AliUEHistograms::SetEtaRange(Float_t etaMin, Float_t etaMax)
1058 {
1059   // sets eta min and max for all contained AliUEHist classes
1060
1061   for (Int_t i=0; i<fgkUEHists; i++)
1062     if (GetUEHist(i))
1063       GetUEHist(i)->SetEtaRange(etaMin, etaMax);
1064 }
1065
1066 //____________________________________________________________________
1067 void AliUEHistograms::SetPtRange(Float_t ptMin, Float_t ptMax)
1068 {
1069   // sets pT min and max for all contained AliUEHist classes
1070
1071   for (Int_t i=0; i<fgkUEHists; i++)
1072     if (GetUEHist(i))
1073       GetUEHist(i)->SetPtRange(ptMin, ptMax);
1074 }
1075
1076 //____________________________________________________________________
1077 void AliUEHistograms::SetPartSpecies(Int_t species)
1078 {
1079   // sets PartSpecie for all contained AliUEHist classes
1080
1081   for (Int_t i=0; i<fgkUEHists; i++)
1082     if (GetUEHist(i))
1083       GetUEHist(i)->SetPartSpecies(species);
1084 }
1085
1086 //____________________________________________________________________
1087 void AliUEHistograms::SetZVtxRange(Float_t min, Float_t max)
1088 {
1089   // sets pT min and max for all contained AliUEHist classes
1090
1091   for (Int_t i=0; i<fgkUEHists; i++)
1092     if (GetUEHist(i))
1093       GetUEHist(i)->SetZVtxRange(min, max);
1094 }
1095
1096 //____________________________________________________________________
1097 void AliUEHistograms::SetContaminationEnhancement(TH1F* hist)
1098 {
1099   // sets the contamination enhancement histogram in all contained AliUEHist classes
1100
1101   for (Int_t i=0; i<fgkUEHists; i++)
1102     if (GetUEHist(i))
1103       GetUEHist(i)->SetContaminationEnhancement(hist);
1104 }
1105
1106 //____________________________________________________________________
1107 void AliUEHistograms::SetCombineMinMax(Bool_t flag)
1108 {
1109   // sets pT min and max for all contained AliUEHist classes
1110
1111   for (Int_t i=0; i<fgkUEHists; i++)
1112     if (GetUEHist(i))
1113       GetUEHist(i)->SetCombineMinMax(flag);
1114 }
1115
1116 //____________________________________________________________________
1117 void AliUEHistograms::SetTrackEtaCut(Float_t value)
1118 {
1119   // sets track eta cut for all contained AliUEHist classes
1120
1121   for (Int_t i=0; i<fgkUEHists; i++)
1122     if (GetUEHist(i))
1123       GetUEHist(i)->SetTrackEtaCut(value);
1124 }
1125
1126 //____________________________________________________________________
1127 void AliUEHistograms::SetWeightPerEvent(Bool_t flag)
1128 {
1129   // sets fWeightPerEvent for all contained AliUEHist classes
1130
1131   fWeightPerEvent = flag;
1132
1133   for (Int_t i=0; i<fgkUEHists; i++)
1134     if (GetUEHist(i))
1135       GetUEHist(i)->SetWeightPerEvent(fWeightPerEvent);
1136 }
1137
1138 //____________________________________________________________________
1139 void AliUEHistograms::Correct(AliUEHistograms* corrections)
1140 {
1141   // corrects the contained histograms by calling AliUEHist::Correct
1142
1143   for (Int_t i=0; i<fgkUEHists; i++)
1144     if (GetUEHist(i))
1145       GetUEHist(i)->Correct(corrections->GetUEHist(i));
1146 }
1147
1148 //____________________________________________________________________
1149 AliUEHistograms &AliUEHistograms::operator=(const AliUEHistograms &c)
1150 {
1151   // assigment operator
1152
1153   DeleteContainers();
1154
1155   if (this != &c)
1156     ((AliUEHistograms &) c).Copy(*this);
1157
1158   return *this;
1159 }
1160
1161 //____________________________________________________________________
1162 void AliUEHistograms::Copy(TObject& c) const
1163 {
1164   // copy function
1165
1166   AliUEHistograms& target = (AliUEHistograms &) c;
1167
1168   if (fNumberDensitypT)
1169     target.fNumberDensitypT = dynamic_cast<AliUEHist*> (fNumberDensitypT->Clone());
1170
1171   if (fSumpT)
1172     target.fSumpT = dynamic_cast<AliUEHist*> (fSumpT->Clone());
1173
1174   if (fNumberDensityPhi)
1175     target.fNumberDensityPhi = dynamic_cast<AliUEHist*> (fNumberDensityPhi->Clone());
1176
1177   if (fCorrelationpT)
1178     target.fCorrelationpT = dynamic_cast<TH2F*> (fCorrelationpT->Clone());
1179
1180   if (fCorrelationEta)
1181     target.fCorrelationEta = dynamic_cast<TH2F*> (fCorrelationEta->Clone());
1182
1183   if (fCorrelationPhi)
1184     target.fCorrelationPhi = dynamic_cast<TH2F*> (fCorrelationPhi->Clone());
1185
1186   if (fCorrelationR)
1187     target.fCorrelationR = dynamic_cast<TH2F*> (fCorrelationR->Clone());
1188
1189   if (fCorrelationLeading2Phi)
1190     target.fCorrelationLeading2Phi = dynamic_cast<TH2F*> (fCorrelationLeading2Phi->Clone());
1191
1192   if (fCorrelationMultiplicity)
1193     target.fCorrelationMultiplicity = dynamic_cast<TH2F*> (fCorrelationMultiplicity->Clone());
1194
1195   if (fYields)
1196     target.fYields = dynamic_cast<TH3F*> (fYields->Clone());
1197
1198   if (fEventCount)
1199     target.fEventCount = dynamic_cast<TH2F*> (fEventCount->Clone());
1200
1201   if (fEventCountDifferential)
1202     target.fEventCountDifferential = dynamic_cast<TH3F*> (fEventCountDifferential->Clone());
1203
1204   if (fVertexContributors)
1205     target.fVertexContributors = dynamic_cast<TH1F*> (fVertexContributors->Clone());
1206
1207   if (fCentralityDistribution)
1208     target.fCentralityDistribution = dynamic_cast<TH1F*> (fCentralityDistribution->Clone());
1209
1210   if (fCentralityCorrelation)
1211     target.fCentralityCorrelation = dynamic_cast<TH2F*> (fCentralityCorrelation->Clone());
1212
1213   if (fITSClusterMap)
1214     target.fITSClusterMap = dynamic_cast<TH3F*> (fITSClusterMap->Clone());
1215
1216   if (fControlConvResoncances)
1217     target.fControlConvResoncances = dynamic_cast<TH2F*> (fControlConvResoncances->Clone());
1218
1219   for (Int_t i=0; i<2; i++)
1220     if (fTwoTrackDistancePt[i])
1221       target.fTwoTrackDistancePt[i] = dynamic_cast<TH3F*> (fTwoTrackDistancePt[i]->Clone());
1222
1223   if (fEfficiencyCorrectionTriggers)
1224     target.fEfficiencyCorrectionTriggers = dynamic_cast<THnF*> (fEfficiencyCorrectionTriggers->Clone());
1225
1226  if (fEfficiencyCorrectionAssociated)
1227     target.fEfficiencyCorrectionAssociated = dynamic_cast<THnF*> (fEfficiencyCorrectionAssociated->Clone());
1228
1229   target.fSelectCharge = fSelectCharge;
1230   target.fTriggerSelectCharge = fTriggerSelectCharge;
1231   target.fAssociatedSelectCharge = fAssociatedSelectCharge;
1232   target.fTriggerRestrictEta = fTriggerRestrictEta;
1233   target.fEtaOrdering = fEtaOrdering;
1234   target.fCutConversions = fCutConversions;
1235   target.fCutResonances = fCutResonances;
1236   target.fOnlyOneEtaSide = fOnlyOneEtaSide;
1237   target.fWeightPerEvent = fWeightPerEvent;
1238   target.fRunNumber = fRunNumber;
1239   target.fMergeCount = fMergeCount;
1240   target.fWeightPerEvent = fWeightPerEvent;
1241   target.fPtOrder = fPtOrder;
1242   target.fTwoTrackCutMinRadius = fTwoTrackCutMinRadius;
1243 }
1244
1245 //____________________________________________________________________
1246 Long64_t AliUEHistograms::Merge(TCollection* list)
1247 {
1248   // Merge a list of AliUEHistograms objects with this (needed for
1249   // PROOF).
1250   // Returns the number of merged objects (including this).
1251
1252   if (!list)
1253     return 0;
1254
1255   if (list->IsEmpty())
1256     return 1;
1257
1258   TIterator* iter = list->MakeIterator();
1259   TObject* obj;
1260
1261   // collections of objects
1262   const Int_t kMaxLists = 19;
1263   TList* lists[kMaxLists];
1264
1265   for (Int_t i=0; i<kMaxLists; i++)
1266     lists[i] = new TList;
1267
1268   Int_t count = 0;
1269   while ((obj = iter->Next())) {
1270
1271     AliUEHistograms* entry = dynamic_cast<AliUEHistograms*> (obj);
1272     if (entry == 0)
1273       continue;
1274
1275     if (entry->fNumberDensitypT)
1276       lists[0]->Add(entry->fNumberDensitypT);
1277     if (entry->fSumpT)
1278       lists[1]->Add(entry->fSumpT);
1279     if (entry->fNumberDensityPhi)
1280       lists[2]->Add(entry->fNumberDensityPhi);
1281     lists[3]->Add(entry->fCorrelationpT);
1282     lists[4]->Add(entry->fCorrelationEta);
1283     lists[5]->Add(entry->fCorrelationPhi);
1284     lists[6]->Add(entry->fCorrelationR);
1285     lists[7]->Add(entry->fCorrelationLeading2Phi);
1286     lists[8]->Add(entry->fCorrelationMultiplicity);
1287     lists[9]->Add(entry->fEventCount);
1288     lists[10]->Add(entry->fEventCountDifferential);
1289     lists[11]->Add(entry->fVertexContributors);
1290     lists[12]->Add(entry->fCentralityDistribution);
1291     lists[13]->Add(entry->fITSClusterMap);
1292     if (entry->fTwoTrackDistancePt[0])
1293       lists[14]->Add(entry->fTwoTrackDistancePt[0]);
1294     if (entry->fTwoTrackDistancePt[1])
1295       lists[15]->Add(entry->fTwoTrackDistancePt[1]);
1296     if (entry->fCentralityCorrelation)
1297       lists[16]->Add(entry->fCentralityCorrelation);
1298     if (entry->fYields)
1299       lists[17]->Add(entry->fYields);
1300     if (entry->fControlConvResoncances)
1301       lists[18]->Add(entry->fControlConvResoncances);
1302
1303     fMergeCount += entry->fMergeCount;
1304
1305     count++;
1306   }
1307
1308   if (fNumberDensitypT)
1309     fNumberDensitypT->Merge(lists[0]);
1310   if (fSumpT)
1311     fSumpT->Merge(lists[1]);
1312   if (fNumberDensityPhi)
1313     fNumberDensityPhi->Merge(lists[2]);
1314   fCorrelationpT->Merge(lists[3]);
1315   fCorrelationEta->Merge(lists[4]);
1316   fCorrelationPhi->Merge(lists[5]);
1317   fCorrelationR->Merge(lists[6]);
1318   fCorrelationLeading2Phi->Merge(lists[7]);
1319   fCorrelationMultiplicity->Merge(lists[8]);
1320   fEventCount->Merge(lists[9]);
1321   fEventCountDifferential->Merge(lists[10]);
1322   fVertexContributors->Merge(lists[11]);
1323   fCentralityDistribution->Merge(lists[12]);
1324   fITSClusterMap->Merge(lists[13]);
1325   if (fTwoTrackDistancePt[0] && lists[14]->GetEntries() > 0)
1326     fTwoTrackDistancePt[0]->Merge(lists[14]);
1327   if (fTwoTrackDistancePt[1] && lists[15]->GetEntries() > 0)
1328     fTwoTrackDistancePt[1]->Merge(lists[15]);
1329   if (fCentralityCorrelation)
1330     fCentralityCorrelation->Merge(lists[16]);
1331   if (fYields && lists[17]->GetEntries() > 0)
1332     fYields->Merge(lists[17]);
1333   if (fControlConvResoncances && lists[18]->GetEntries() > 0)
1334     fControlConvResoncances->Merge(lists[18]);
1335
1336   for (Int_t i=0; i<kMaxLists; i++)
1337     delete lists[i];
1338
1339   return count+1;
1340 }
1341
1342 void AliUEHistograms::CopyReconstructedData(AliUEHistograms* from)
1343 {
1344   // copies those histograms extracted from ESD to this object
1345
1346   for (Int_t i=0; i<fgkUEHists; i++)
1347     if (GetUEHist(i))
1348       GetUEHist(i)->CopyReconstructedData(from->GetUEHist(i));
1349 }
1350
1351 void AliUEHistograms::DeepCopy(AliUEHistograms* from)
1352 {
1353   // copies the entries of this object's members from the object <from> to this object
1354
1355   for (Int_t i=0; i<fgkUEHists; i++)
1356     if (GetUEHist(i) && from->GetUEHist(i))
1357       GetUEHist(i)->DeepCopy(from->GetUEHist(i));
1358 }
1359
1360 void AliUEHistograms::ExtendTrackingEfficiency(Bool_t verbose)
1361 {
1362   // delegates to AliUEHists
1363
1364   for (Int_t i=0; i<fgkUEHists; i++)
1365     if (GetUEHist(i))
1366       GetUEHist(i)->ExtendTrackingEfficiency(verbose);
1367 }
1368
1369 void AliUEHistograms::Scale(Double_t factor)
1370 {
1371   // scales all contained histograms by the given factor
1372
1373   for (Int_t i=0; i<fgkUEHists; i++)
1374     if (GetUEHist(i))
1375       GetUEHist(i)->Scale(factor);
1376
1377   TList list;
1378   list.Add(fCorrelationpT);
1379   list.Add(fCorrelationEta);
1380   list.Add(fCorrelationPhi);
1381   list.Add(fCorrelationR);
1382   list.Add(fCorrelationLeading2Phi);
1383   list.Add(fCorrelationMultiplicity);
1384   list.Add(fYields);
1385   list.Add(fEventCount);
1386   list.Add(fEventCountDifferential);
1387   list.Add(fVertexContributors);
1388   list.Add(fCentralityDistribution);
1389   list.Add(fCentralityCorrelation);
1390   list.Add(fITSClusterMap);
1391   list.Add(fTwoTrackDistancePt[0]);
1392   list.Add(fTwoTrackDistancePt[1]);
1393   list.Add(fControlConvResoncances);
1394
1395   for (Int_t i=0; i<list.GetEntries(); i++)
1396     ((TH1*) list.At(i))->Scale(factor);
1397 }
1398
1399 void AliUEHistograms::Reset()
1400 {
1401   // delegates to AliUEHists
1402
1403   for (Int_t i=0; i<fgkUEHists; i++)
1404     if (GetUEHist(i))
1405       GetUEHist(i)->Reset();
1406 }