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 "TH2F.h"
  33 #include "TH1F.h"
  34 #include "TH3F.h"
  35 #include "TMath.h"
  36 #include "TLorentzVector.h"
  37
  38 ClassImp(AliUEHistograms)
  39
  40 const Int_t AliUEHistograms::fgkUEHists = 3;
  41
  42 AliUEHistograms::AliUEHistograms(const char* name, const char* histograms) :
  43   TNamed(name, name),
  44   fNumberDensitypT(0),
  45   fSumpT(0),
  46   fNumberDensityPhi(0),
  47   fCorrelationpT(0),
  48   fCorrelationEta(0),
  49   fCorrelationPhi(0),
  50   fCorrelationR(0),
  51   fCorrelationLeading2Phi(0),
  52   fCorrelationMultiplicity(0),
  53   fYields(0),
  54   fEventCount(0),
  55   fEventCountDifferential(0),
  56   fVertexContributors(0),
  57   fCentralityDistribution(0),
  58   fCentralityCorrelation(0),
  59   fITSClusterMap(0),
  60   fSelectCharge(0),
  61   fTriggerRestrictEta(-1),
  62   fEtaOrdering(kFALSE),
  63   fCutConversions(kFALSE),
  64   fCutResonances(kFALSE),
  65   fOnlyOneEtaSide(0),
  66   fRunNumber(0)
  67 {
  68   // Constructor
  69   //
  70   // the string histograms defines which histograms are created:
  71   //    1 = NumberDensitypT
  72   //    2 = SumpT
  73   //    3 = NumberDensityPhi
  74   //    4 = NumberDensityPhiCentrality (other multiplicity for Pb)
  75
  76   AliLog::SetClassDebugLevel("AliCFContainer", -1);
  77   AliLog::SetClassDebugLevel("AliCFGridSparse", -3);
  78
  79   fTwoTrackDistancePt[0] = 0;
  80   fTwoTrackDistancePt[1] = 0;
  81
  82   TString histogramsStr(histograms);
  83
  84   if (histogramsStr.Contains("1"))
  85     fNumberDensitypT = new AliUEHist("NumberDensitypT");
  86   if (histogramsStr.Contains("2"))
  87     fSumpT = new AliUEHist("SumpT");
  88
  89   if (histogramsStr.Contains("3"))
  90     fNumberDensityPhi = new AliUEHist("NumberDensityPhi");
  91   else if (histogramsStr.Contains("4R"))
  92     fNumberDensityPhi = new AliUEHist("NumberDensityPhiCentrality");
  93   else if (histogramsStr.Contains("4"))
  94     fNumberDensityPhi = new AliUEHist("NumberDensityPhiCentralityTTR");
  95   else if (histogramsStr.Contains("5R"))
  96     fNumberDensityPhi = new AliUEHist("NumberDensityPhiCentralityVtx");
  97   else if (histogramsStr.Contains("5"))
  98     fNumberDensityPhi = new AliUEHist("NumberDensityPhiCentralityVtxTTR");
  99
 100   // do not add this hists to the directory
 101   Bool_t oldStatus = TH1::AddDirectoryStatus();
 102   TH1::AddDirectory(kFALSE);
 103
 104   if (!histogramsStr.Contains("4") && !histogramsStr.Contains("5"))
 105   {
 106     fCorrelationpT  = new TH2F("fCorrelationpT", ";p_{T,lead} (MC);p_{T,lead} (RECO)", 200, 0, 50, 200, 0, 50);
 107     fCorrelationEta = new TH2F("fCorrelationEta", ";#eta_{lead} (MC);#eta_{T,lead} (RECO)", 200, -1, 1, 200, -1, 1);
 108     fCorrelationPhi = new TH2F("fCorrelationPhi", ";#varphi_{lead} (MC);#varphi_{T,lead} (RECO)", 200, 0, TMath::TwoPi(), 200, 0, TMath::TwoPi());
 109   }
 110   else
 111   {
 112     fCorrelationpT  = new TH2F("fCorrelationpT", ";Centrality;p_{T} (RECO)", 100, 0, 100.001, 200, 0, 50);
 113     fCorrelationEta = new TH2F("fCorrelationEta", ";Centrality;#eta (RECO)", 100, 0, 100.001, 200, -5, 5);
 114     fCorrelationPhi = new TH2F("fCorrelationPhi", ";Centrality;#varphi (RECO)", 100, 0, 100.001, 200, 0, TMath::TwoPi());
 115   }
 116
 117   fCorrelationR =   new TH2F("fCorrelationR", ";R;p_{T,lead} (MC)", 200, 0, 2, 200, 0, 50);
 118   fCorrelationLeading2Phi = new TH2F("fCorrelationLeading2Phi", ";#Delta #varphi;p_{T,lead} (MC)", 200, -TMath::Pi(), TMath::Pi(), 200, 0, 50);
 119   fCorrelationMultiplicity = new TH2F("fCorrelationMultiplicity", ";MC tracks;Reco tracks", 100, -0.5, 99.5, 100, -0.5, 99.5);
 120   fYields = new TH3F("fYields", ";centrality;pT;eta", 100, 0, 100, 40, 0, 20, 20, -1, 1);
 121
 122   if (!histogramsStr.Contains("4") && !histogramsStr.Contains("5"))
 123   {
 124     fEventCount = new TH2F("fEventCount", ";step;event type;count", AliUEHist::fgkCFSteps+2, -2.5, -0.5 + AliUEHist::fgkCFSteps, 3, -0.5, 2.5);
 125     fEventCount->GetYaxis()->SetBinLabel(1, "ND");
 126     fEventCount->GetYaxis()->SetBinLabel(2, "SD");
 127     fEventCount->GetYaxis()->SetBinLabel(3, "DD");
 128   }
 129   else
 130   {
 131     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());
 132   }
 133
 134   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);
 135   fEventCountDifferential->GetZaxis()->SetBinLabel(1, "ND");
 136   fEventCountDifferential->GetZaxis()->SetBinLabel(2, "SD");
 137   fEventCountDifferential->GetZaxis()->SetBinLabel(3, "DD");
 138
 139   fVertexContributors = new TH1F("fVertexContributors", ";contributors;count", 100, -0.5, 99.5);
 140
 141   if (fNumberDensityPhi)
 142   {
 143     fCentralityDistribution = new TH1F("fCentralityDistribution", ";centrality;count", fNumberDensityPhi->GetEventHist()->GetNBins(1), fNumberDensityPhi->GetEventHist()->GetAxis(1, 0)->GetXbins()->GetArray());
 144     fCentralityCorrelation = new TH2F("fCentralityCorrelation", ";centrality;multiplicity", 101, 0, 101, 200, 0, 4000);
 145   }
 146
 147   fITSClusterMap = new TH3F("fITSClusterMap", "; its cluster map; centrality; pT", 256, -0.5, 255.5, 20, 0, 100.001, 100, 0, 20);
 148
 149   TH1::AddDirectory(oldStatus);
 150 }
 151
 152 //_____________________________________________________________________________
 153 AliUEHistograms::AliUEHistograms(const AliUEHistograms &c) :
 154   TNamed(fName, fTitle),
 155   fNumberDensitypT(0),
 156   fSumpT(0),
 157   fNumberDensityPhi(0),
 158   fCorrelationpT(0),
 159   fCorrelationEta(0),
 160   fCorrelationPhi(0),
 161   fCorrelationR(0),
 162   fCorrelationLeading2Phi(0),
 163   fCorrelationMultiplicity(0),
 164   fYields(0),
 165   fEventCount(0),
 166   fEventCountDifferential(0),
 167   fVertexContributors(0),
 168   fCentralityDistribution(0),
 169   fCentralityCorrelation(0),
 170   fITSClusterMap(0),
 171   fSelectCharge(0),
 172   fTriggerRestrictEta(-1),
 173   fEtaOrdering(kFALSE),
 174   fCutConversions(kFALSE),
 175   fCutResonances(kFALSE),
 176   fOnlyOneEtaSide(0),
 177   fRunNumber(0)
 178 {
 179   //
 180   // AliUEHistograms copy constructor
 181   //
 182
 183   fTwoTrackDistancePt[0] = 0;
 184   fTwoTrackDistancePt[1] = 0;
 185
 186   ((AliUEHistograms &) c).Copy(*this);
 187 }
 188
 189 //____________________________________________________________________
 190 AliUEHistograms::~AliUEHistograms()
 191 {
 192   // Destructor
 193
 194   DeleteContainers();
 195 }
 196
 197 void AliUEHistograms::DeleteContainers()
 198 {
 199   if (fNumberDensitypT)
 200   {
 201     delete fNumberDensitypT;
 202     fNumberDensitypT = 0;
 203   }
 204
 205   if (fSumpT)
 206   {
 207     delete fSumpT;
 208     fSumpT = 0;
 209   }
 210
 211   if (fNumberDensityPhi)
 212   {
 213     delete fNumberDensityPhi;
 214     fNumberDensityPhi = 0;
 215   }
 216
 217   if (fCorrelationpT)
 218   {
 219     delete fCorrelationpT;
 220     fCorrelationpT = 0;
 221   }
 222
 223   if (fCorrelationEta)
 224   {
 225     delete fCorrelationEta;
 226     fCorrelationEta = 0;
 227   }
 228
 229   if (fCorrelationPhi)
 230   {
 231     delete fCorrelationPhi;
 232     fCorrelationPhi = 0;
 233   }
 234
 235   if (fCorrelationR)
 236   {
 237     delete fCorrelationR;
 238     fCorrelationR = 0;
 239   }
 240
 241   if (fCorrelationLeading2Phi)
 242   {
 243     delete fCorrelationLeading2Phi;
 244     fCorrelationLeading2Phi = 0;
 245   }
 246
 247   if (fCorrelationMultiplicity)
 248   {
 249     delete fCorrelationMultiplicity;
 250     fCorrelationMultiplicity = 0;
 251   }
 252
 253   if (fYields)
 254   {
 255     delete fYields;
 256     fYields = 0;
 257   }
 258
 259   if (fEventCount)
 260   {
 261     delete fEventCount;
 262     fEventCount = 0;
 263   }
 264
 265   if (fEventCountDifferential)
 266   {
 267     delete fEventCountDifferential;
 268     fEventCountDifferential = 0;
 269   }
 270
 271   if (fVertexContributors)
 272   {
 273     delete fVertexContributors;
 274     fVertexContributors = 0;
 275   }
 276
 277   if (fCentralityDistribution)
 278   {
 279     delete fCentralityDistribution;
 280     fCentralityDistribution = 0;
 281   }
 282
 283   if (fCentralityCorrelation)
 284   {
 285     delete fCentralityCorrelation;
 286     fCentralityCorrelation = 0;
 287   }
 288
 289   if (fITSClusterMap)
 290   {
 291     delete fITSClusterMap;
 292     fITSClusterMap = 0;
 293   }
 294
 295   for (Int_t i=0; i<2; i++)
 296     if (fTwoTrackDistancePt[i])
 297     {
 298       delete fTwoTrackDistancePt[i];
 299       fTwoTrackDistancePt[i] = 0;
 300     }
 301 }
 302
 303 AliUEHist* AliUEHistograms::GetUEHist(Int_t id)
 304 {
 305   // returns AliUEHist object, useful for loops
 306
 307   switch (id)
 308   {
 309     case 0: return fNumberDensitypT; break;
 310     case 1: return fSumpT; break;
 311     case 2: return fNumberDensityPhi; break;
 312   }
 313
 314   return 0;
 315 }
 316
 317 //____________________________________________________________________
 318 Int_t AliUEHistograms::CountParticles(TList* list, Float_t ptMin)
 319 {
 320   // counts the number of particles in the list with a pT above ptMin
 321   // TODO eta cut needed here?
 322
 323   Int_t count = 0;
 324   for (Int_t j=0; j<list->GetEntries(); j++)
 325     if (((AliVParticle*) list->At(j))->Pt() > ptMin)
 326       count++;
 327
 328   return count;
 329 }
 330
 331 //____________________________________________________________________
 332 void AliUEHistograms::Fill(Int_t eventType, Float_t zVtx, AliUEHist::CFStep step, AliVParticle* leading, TList* toward, TList* away, TList* min, TList* max)
 333 {
 334   // fills the UE event histograms
 335   //
 336   // 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
 337
 338   // if leading is not set, just fill event statistics
 339   if (leading)
 340   {
 341     Int_t multiplicity = 0;
 342
 343     // TODO configurable?
 344     Float_t ptMin = 0.15;
 345     if (leading->Pt() > ptMin)
 346       multiplicity++;
 347
 348     multiplicity += CountParticles(toward, ptMin);
 349     multiplicity += CountParticles(away, ptMin);
 350     multiplicity += CountParticles(min, ptMin);
 351     multiplicity += CountParticles(max, ptMin);
 352
 353     FillRegion(AliUEHist::kToward, zVtx, step, leading, toward, multiplicity);
 354     FillRegion(AliUEHist::kAway,   zVtx, step, leading, away, multiplicity);
 355     FillRegion(AliUEHist::kMin,    zVtx, step, leading, min, multiplicity);
 356     FillRegion(AliUEHist::kMax,    zVtx, step, leading, max, multiplicity);
 357
 358     Double_t vars[3];
 359     vars[0] = leading->Pt();
 360     vars[1] = multiplicity;
 361     vars[2] = zVtx;
 362     for (Int_t i=0; i<fgkUEHists; i++)
 363       if (GetUEHist(i))
 364         GetUEHist(i)->GetEventHist()->Fill(vars, step);
 365
 366     fEventCountDifferential->Fill(leading->Pt(), step, eventType);
 367   }
 368
 369   FillEvent(eventType, step);
 370 }
 371
 372 //____________________________________________________________________
 373 void AliUEHistograms::FillRegion(AliUEHist::Region region, Float_t zVtx, AliUEHist::CFStep step, AliVParticle* leading, TList* list, Int_t multiplicity)
 374 {
 375   // loops over AliVParticles in list and fills the given region at the given step
 376   //
 377   // See also Fill(...)
 378
 379   for (Int_t i=0; i<list->GetEntries(); i++)
 380   {
 381     AliVParticle* particle = (AliVParticle*) list->At(i);
 382
 383     Double_t vars[6];
 384     vars[0] = particle->Eta();
 385     vars[1] = particle->Pt();
 386     vars[2] = leading->Pt();
 387     vars[3] = multiplicity;
 388     vars[4] = leading->Phi() - particle->Phi();
 389     if (vars[4] > 1.5 * TMath::Pi())
 390       vars[4] -= TMath::TwoPi();
 391     if (vars[4] < -0.5 * TMath::Pi())
 392       vars[4] += TMath::TwoPi();
 393     vars[5] = zVtx;
 394
 395     if (fNumberDensitypT)
 396       fNumberDensitypT->GetTrackHist(region)->Fill(vars, step);
 397
 398     if (fSumpT)
 399       fSumpT->GetTrackHist(region)->Fill(vars, step, particle->Pt());
 400
 401     // fill all in toward region (is anyway as function of delta phi!)
 402     if (fNumberDensityPhi)
 403       fNumberDensityPhi->GetTrackHist(AliUEHist::kToward)->Fill(vars, step);
 404   }
 405 }
 406
 407 //____________________________________________________________________
 408 void AliUEHistograms::Fill(AliVParticle* leadingMC, AliVParticle* leadingReco)
 409 {
 410   // fills the correlation histograms
 411
 412   if (leadingMC)
 413   {
 414     fCorrelationpT->Fill(leadingMC->Pt(), leadingReco->Pt());
 415     if (leadingMC->Pt() > 0.5)
 416     {
 417       fCorrelationEta->Fill(leadingMC->Eta(), leadingReco->Eta());
 418       fCorrelationPhi->Fill(leadingMC->Phi(), leadingReco->Phi());
 419     }
 420
 421     Float_t phiDiff = leadingMC->Phi() - leadingReco->Phi();
 422     if (phiDiff > TMath::Pi())
 423       phiDiff -= TMath::TwoPi();
 424     if (phiDiff < -TMath::Pi())
 425       phiDiff += TMath::TwoPi();
 426
 427     Float_t etaDiff = leadingMC->Eta() - leadingReco->Eta();
 428
 429     fCorrelationR->Fill(TMath::Sqrt(phiDiff * phiDiff + etaDiff * etaDiff), leadingMC->Pt());
 430     fCorrelationLeading2Phi->Fill(phiDiff, leadingMC->Pt());
 431   }
 432   else
 433   {
 434     fCorrelationpT->Fill(1.0, leadingReco->Pt());
 435     if (leadingReco->Pt() > 0.5)
 436     {
 437       fCorrelationEta->Fill(0.0, leadingReco->Eta());
 438       fCorrelationPhi->Fill(0.0, leadingReco->Phi());
 439     }
 440   }
 441 }
 442
 443 //____________________________________________________________________
 444 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)
 445 {
 446   // fills the fNumberDensityPhi histogram
 447   //
 448   // this function need a list of AliVParticles which contain the particles/tracks to be filled
 449   //
 450   // if mixed is non-0, mixed events are filled, the trigger particle is from particles, the associated from mixed
 451   // if weight < 0, then the pt of the associated particle is filled as weight
 452
 453   Bool_t fillpT = kFALSE;
 454   if (weight < 0)
 455     fillpT = kTRUE;
 456
 457   if (twoTrackEfficiencyCut && !fTwoTrackDistancePt[0])
 458   {
 459     // do not add this hists to the directory
 460     Bool_t oldStatus = TH1::AddDirectoryStatus();
 461     TH1::AddDirectory(kFALSE);
 462
 463     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);
 464     fTwoTrackDistancePt[1] = (TH3F*) fTwoTrackDistancePt[0]->Clone("fTwoTrackDistancePt[1]");
 465
 466     TH1::AddDirectory(oldStatus);
 467   }
 468
 469   // Eta() is extremely time consuming, therefore cache it for the inner loop here:
 470   TObjArray* input = (mixed) ? mixed : particles;
 471   TArrayF eta(input->GetEntriesFast());
 472   for (Int_t i=0; i<input->GetEntriesFast(); i++)
 473     eta[i] = ((AliVParticle*) input->At(i))->Eta();
 474
 475   // if particles is not set, just fill event statistics
 476   if (particles)
 477   {
 478     Int_t jMax = particles->GetEntriesFast();
 479     if (mixed)
 480       jMax = mixed->GetEntriesFast();
 481
 482     for (Int_t i=0; i<particles->GetEntriesFast(); i++)
 483     {
 484       AliVParticle* triggerParticle = (AliVParticle*) particles->At(i);
 485
 486       // some optimization
 487       Float_t triggerEta = triggerParticle->Eta();
 488
 489       if (fTriggerRestrictEta > 0 && TMath::Abs(triggerEta) > fTriggerRestrictEta)
 490         continue;
 491
 492       if (fOnlyOneEtaSide != 0)
 493       {
 494         if (fOnlyOneEtaSide * triggerEta < 0)
 495           continue;
 496       }
 497
 498       if (!mixed)
 499       {
 500         // QA
 501         fCorrelationpT->Fill(centrality, triggerParticle->Pt());
 502         fCorrelationEta->Fill(centrality, triggerEta);
 503         fCorrelationPhi->Fill(centrality, triggerParticle->Phi());
 504         fYields->Fill(centrality, triggerParticle->Pt(), triggerEta);
 505 /*        if (dynamic_cast<AliAODTrack*>(triggerParticle))
 506           fITSClusterMap->Fill(((AliAODTrack*) triggerParticle)->GetITSClusterMap(), centrality, triggerParticle->Pt());*/
 507       }
 508
 509       for (Int_t j=0; j<jMax; j++)
 510       {
 511         if (!mixed && i == j)
 512           continue;
 513
 514         AliVParticle* particle = 0;
 515         if (!mixed)
 516           particle = (AliVParticle*) particles->At(j);
 517         else
 518           particle = (AliVParticle*) mixed->At(j);
 519
 520         // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event for cross-checks)
 521         if (mixed && triggerParticle == particle)
 522           continue;
 523
 524         if (particle->Pt() > triggerParticle->Pt())
 525           continue;
 526
 527         if (fSelectCharge > 0)
 528         {
 529           // skip like sign
 530           if (fSelectCharge == 1 && particle->Charge() * triggerParticle->Charge() > 0)
 531             continue;
 532
 533           // skip unlike sign
 534           if (fSelectCharge == 2 && particle->Charge() * triggerParticle->Charge() < 0)
 535             continue;
 536         }
 537
 538         if (fEtaOrdering)
 539         {
 540           if (triggerEta < 0 && eta[j] < triggerEta)
 541             continue;
 542           if (triggerEta > 0 && eta[j] > triggerEta)
 543             continue;
 544         }
 545
 546         // conversions
 547         if (fCutConversions && particle->Charge() * triggerParticle->Charge() < 0)
 548         {
 549           Float_t mass = GetInvMassSquared(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), particle->Pt(), eta[j], particle->Phi(), 0.510e-3);
 550
 551           if (mass < 0.04*0.04)
 552             continue;
 553         }
 554
 555         // K0s, rhos
 556         if (fCutResonances && particle->Charge() * triggerParticle->Charge() < 0)
 557         {
 558           Float_t mass = GetInvMassSquared(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), particle->Pt(), eta[j], particle->Phi(), 0.1396);
 559
 560           if ((mass > 0.49*0.49 && mass < 0.51*0.51) || (mass > 0.765*0.765 && mass < 0.785*0.785))
 561             continue;
 562         }
 563
 564         if (twoTrackEfficiencyCut)
 565         {
 566           // the variables & cuthave been developed by the HBT group
 567           // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
 568
 569           Float_t phi1 = triggerParticle->Phi();
 570           Float_t pt1 = triggerParticle->Pt();
 571           Float_t charge1 = triggerParticle->Charge();
 572
 573           Float_t phi2 = particle->Phi();
 574           Float_t pt2 = particle->Pt();
 575           Float_t charge2 = particle->Charge();
 576
 577           Float_t deta = triggerEta - eta[j];
 578
 579           // optimization
 580           if (TMath::Abs(deta) < twoTrackEfficiencyCutValue * 2.5 * 3)
 581           {
 582             // check first boundaries to see if is worth to loop and find the minimum
 583             Float_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 0.8, bSign);
 584             Float_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 2.5, bSign);
 585
 586             const Float_t kLimit = twoTrackEfficiencyCutValue * 3;
 587
 588             Float_t dphistarminabs = 1e5;
 589             Float_t dphistarmin = 1e5;
 590             if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
 591             {
 592               for (Double_t rad=0.8; rad<2.51; rad+=0.01)
 593               {
 594                 Float_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
 595
 596                 Float_t dphistarabs = TMath::Abs(dphistar);
 597
 598                 if (dphistarabs < dphistarminabs)
 599                 {
 600                   dphistarmin = dphistar;
 601                   dphistarminabs = dphistarabs;
 602                 }
 603               }
 604
 605               fTwoTrackDistancePt[0]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));
 606
 607               if (dphistarminabs < twoTrackEfficiencyCutValue && TMath::Abs(deta) < twoTrackEfficiencyCutValue)
 608               {
 609 //              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);
 610                 continue;
 611               }
 612
 613               fTwoTrackDistancePt[1]->Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2));
 614             }
 615           }
 616         }
 617
 618         Double_t vars[6];
 619         vars[0] = triggerEta - eta[j];
 620         vars[1] = particle->Pt();
 621         vars[2] = triggerParticle->Pt();
 622         vars[3] = centrality;
 623         vars[4] = triggerParticle->Phi() - particle->Phi();
 624         if (vars[4] > 1.5 * TMath::Pi())
 625           vars[4] -= TMath::TwoPi();
 626         if (vars[4] < -0.5 * TMath::Pi())
 627           vars[4] += TMath::TwoPi();
 628         vars[5] = zVtx;
 629
 630         if (fillpT)
 631           weight = particle->Pt();
 632
 633         // fill all in toward region and do not use the other regions
 634         fNumberDensityPhi->GetTrackHist(AliUEHist::kToward)->Fill(vars, step, weight);
 635
 636 //      Printf("%.2f %.2f --> %.2f", triggerEta, eta[j], vars[0]);
 637       }
 638
 639       if (firstTime)
 640       {
 641         // once per trigger particle
 642         Double_t vars[3];
 643         vars[0] = triggerParticle->Pt();
 644         vars[1] = centrality;
 645         vars[2] = zVtx;
 646         fNumberDensityPhi->GetEventHist()->Fill(vars, step);
 647       }
 648     }
 649   }
 650
 651   fCentralityDistribution->Fill(centrality);
 652   fCentralityCorrelation->Fill(centrality, particles->GetEntriesFast());
 653   FillEvent(centrality, step);
 654 }
 655
 656 //____________________________________________________________________
 657 void AliUEHistograms::FillTrackingEfficiency(TObjArray* mc, TObjArray* recoPrim, TObjArray* recoAll, Int_t particleType, Double_t centrality)
 658 {
 659   // fills the tracking efficiency objects
 660   //
 661   // mc: all primary MC particles
 662   // recoPrim: reconstructed primaries (again MC particles)
 663   // recoAll: reconstructed (again MC particles)
 664   // particleType is: 0 for pion, 1 for kaon, 2 for proton, 3 for others
 665
 666   for (Int_t step=0; step<3; step++)
 667   {
 668     TObjArray* list = mc;
 669     if (step == 1)
 670       list = recoPrim;
 671     else if (step == 2)
 672       list = recoAll;
 673
 674     for (Int_t i=0; i<list->GetEntriesFast(); i++)
 675     {
 676       AliVParticle* particle = (AliVParticle*) list->At(i);
 677       Double_t vars[4];
 678       vars[0] = particle->Eta();
 679       vars[1] = particle->Pt();
 680       vars[2] = particleType;
 681       vars[3] = centrality;
 682
 683       for (Int_t j=0; j<fgkUEHists; j++)
 684         if (GetUEHist(j))
 685           GetUEHist(j)->GetTrackHistEfficiency()->Fill(vars, step);
 686     }
 687   }
 688 }
 689
 690 //____________________________________________________________________
 691 void AliUEHistograms::FillEvent(Int_t eventType, Int_t step)
 692 {
 693   // fills the number of events at the given step and the given enty type
 694   //
 695   // WARNING: This function is called from Fill, so only call it for steps where Fill is not called
 696
 697   fEventCount->Fill(step, eventType);
 698 }
 699
 700 //____________________________________________________________________
 701 void AliUEHistograms::FillEvent(Double_t centrality, Int_t step)
 702 {
 703   // fills the number of events at the given step and the given centrality
 704   //
 705   // WARNING: This function is called from Fill, so only call it for steps where Fill is not called
 706
 707   fEventCount->Fill(step, centrality);
 708 }
 709
 710 //____________________________________________________________________
 711 void AliUEHistograms::SetEtaRange(Float_t etaMin, Float_t etaMax)
 712 {
 713   // sets eta min and max for all contained AliUEHist classes
 714
 715   for (Int_t i=0; i<fgkUEHists; i++)
 716     if (GetUEHist(i))
 717       GetUEHist(i)->SetEtaRange(etaMin, etaMax);
 718 }
 719
 720 //____________________________________________________________________
 721 void AliUEHistograms::SetPtRange(Float_t ptMin, Float_t ptMax)
 722 {
 723   // sets pT min and max for all contained AliUEHist classes
 724
 725   for (Int_t i=0; i<fgkUEHists; i++)
 726     if (GetUEHist(i))
 727       GetUEHist(i)->SetPtRange(ptMin, ptMax);
 728 }
 729
 730 //____________________________________________________________________
 731 void AliUEHistograms::SetZVtxRange(Float_t min, Float_t max)
 732 {
 733   // sets pT min and max for all contained AliUEHist classes
 734
 735   for (Int_t i=0; i<fgkUEHists; i++)
 736     if (GetUEHist(i))
 737       GetUEHist(i)->SetZVtxRange(min, max);
 738 }
 739
 740 //____________________________________________________________________
 741 void AliUEHistograms::SetContaminationEnhancement(TH1F* hist)
 742 {
 743   // sets the contamination enhancement histogram in all contained AliUEHist classes
 744
 745   for (Int_t i=0; i<fgkUEHists; i++)
 746     if (GetUEHist(i))
 747       GetUEHist(i)->SetContaminationEnhancement(hist);
 748 }
 749
 750 //____________________________________________________________________
 751 void AliUEHistograms::SetCombineMinMax(Bool_t flag)
 752 {
 753   // sets pT min and max for all contained AliUEHist classes
 754
 755   for (Int_t i=0; i<fgkUEHists; i++)
 756     if (GetUEHist(i))
 757       GetUEHist(i)->SetCombineMinMax(flag);
 758 }
 759
 760 //____________________________________________________________________
 761 void AliUEHistograms::Correct(AliUEHistograms* corrections)
 762 {
 763   // corrects the contained histograms by calling AliUEHist::Correct
 764
 765   for (Int_t i=0; i<fgkUEHists; i++)
 766     if (GetUEHist(i))
 767       GetUEHist(i)->Correct(corrections->GetUEHist(i));
 768 }
 769
 770 //____________________________________________________________________
 771 AliUEHistograms &AliUEHistograms::operator=(const AliUEHistograms &c)
 772 {
 773   // assigment operator
 774
 775   DeleteContainers();
 776
 777   if (this != &c)
 778     ((AliUEHistograms &) c).Copy(*this);
 779
 780   return *this;
 781 }
 782
 783 //____________________________________________________________________
 784 void AliUEHistograms::Copy(TObject& c) const
 785 {
 786   // copy function
 787
 788   AliUEHistograms& target = (AliUEHistograms &) c;
 789
 790   if (fNumberDensitypT)
 791     target.fNumberDensitypT = dynamic_cast<AliUEHist*> (fNumberDensitypT->Clone());
 792
 793   if (fSumpT)
 794     target.fSumpT = dynamic_cast<AliUEHist*> (fSumpT->Clone());
 795
 796   if (fNumberDensityPhi)
 797     target.fNumberDensityPhi = dynamic_cast<AliUEHist*> (fNumberDensityPhi->Clone());
 798
 799   if (fCorrelationpT)
 800     target.fCorrelationpT = dynamic_cast<TH2F*> (fCorrelationpT->Clone());
 801
 802   if (fCorrelationEta)
 803     target.fCorrelationEta = dynamic_cast<TH2F*> (fCorrelationEta->Clone());
 804
 805   if (fCorrelationPhi)
 806     target.fCorrelationPhi = dynamic_cast<TH2F*> (fCorrelationPhi->Clone());
 807
 808   if (fCorrelationR)
 809     target.fCorrelationR = dynamic_cast<TH2F*> (fCorrelationR->Clone());
 810
 811   if (fCorrelationLeading2Phi)
 812     target.fCorrelationLeading2Phi = dynamic_cast<TH2F*> (fCorrelationLeading2Phi->Clone());
 813
 814   if (fCorrelationMultiplicity)
 815     target.fCorrelationMultiplicity = dynamic_cast<TH2F*> (fCorrelationMultiplicity->Clone());
 816
 817   if (fYields)
 818     target.fYields = dynamic_cast<TH3F*> (fYields->Clone());
 819
 820   if (fEventCount)
 821     target.fEventCount = dynamic_cast<TH2F*> (fEventCount->Clone());
 822
 823   if (fEventCountDifferential)
 824     target.fEventCountDifferential = dynamic_cast<TH3F*> (fEventCountDifferential->Clone());
 825
 826   if (fVertexContributors)
 827     target.fVertexContributors = dynamic_cast<TH1F*> (fVertexContributors->Clone());
 828
 829   if (fCentralityDistribution)
 830     target.fCentralityDistribution = dynamic_cast<TH1F*> (fCentralityDistribution->Clone());
 831
 832   if (fCentralityCorrelation)
 833     target.fCentralityCorrelation = dynamic_cast<TH2F*> (fCentralityCorrelation->Clone());
 834
 835   if (fITSClusterMap)
 836     target.fITSClusterMap = dynamic_cast<TH3F*> (fITSClusterMap->Clone());
 837
 838   for (Int_t i=0; i<2; i++)
 839     if (fTwoTrackDistancePt[i])
 840       target.fTwoTrackDistancePt[i] = dynamic_cast<TH3F*> (fTwoTrackDistancePt[i]->Clone());
 841
 842   target.fSelectCharge = fSelectCharge;
 843   target.fTriggerRestrictEta = fTriggerRestrictEta;
 844   target.fEtaOrdering = fEtaOrdering;
 845   target.fCutConversions = fCutConversions;
 846   target.fCutResonances = fCutResonances;
 847   target.fOnlyOneEtaSide = fOnlyOneEtaSide;
 848   target.fRunNumber = fRunNumber;
 849 }
 850
 851 //____________________________________________________________________
 852 Long64_t AliUEHistograms::Merge(TCollection* list)
 853 {
 854   // Merge a list of AliUEHistograms objects with this (needed for
 855   // PROOF).
 856   // Returns the number of merged objects (including this).
 857
 858   if (!list)
 859     return 0;
 860
 861   if (list->IsEmpty())
 862     return 1;
 863
 864   TIterator* iter = list->MakeIterator();
 865   TObject* obj;
 866
 867   // collections of objects
 868   const Int_t kMaxLists = 18;
 869   TList* lists[kMaxLists];
 870
 871   for (Int_t i=0; i<kMaxLists; i++)
 872     lists[i] = new TList;
 873
 874   Int_t count = 0;
 875   while ((obj = iter->Next())) {
 876
 877     AliUEHistograms* entry = dynamic_cast<AliUEHistograms*> (obj);
 878     if (entry == 0)
 879       continue;
 880
 881     if (entry->fNumberDensitypT)
 882       lists[0]->Add(entry->fNumberDensitypT);
 883     if (entry->fSumpT)
 884       lists[1]->Add(entry->fSumpT);
 885     if (entry->fNumberDensityPhi)
 886       lists[2]->Add(entry->fNumberDensityPhi);
 887     lists[3]->Add(entry->fCorrelationpT);
 888     lists[4]->Add(entry->fCorrelationEta);
 889     lists[5]->Add(entry->fCorrelationPhi);
 890     lists[6]->Add(entry->fCorrelationR);
 891     lists[7]->Add(entry->fCorrelationLeading2Phi);
 892     lists[8]->Add(entry->fCorrelationMultiplicity);
 893     lists[9]->Add(entry->fEventCount);
 894     lists[10]->Add(entry->fEventCountDifferential);
 895     lists[11]->Add(entry->fVertexContributors);
 896     lists[12]->Add(entry->fCentralityDistribution);
 897     lists[13]->Add(entry->fITSClusterMap);
 898     if (entry->fTwoTrackDistancePt[0])
 899       lists[14]->Add(entry->fTwoTrackDistancePt[0]);
 900     if (entry->fTwoTrackDistancePt[1])
 901       lists[15]->Add(entry->fTwoTrackDistancePt[1]);
 902     if (entry->fCentralityCorrelation)
 903       lists[16]->Add(entry->fCentralityCorrelation);
 904     if (entry->fYields)
 905       lists[17]->Add(entry->fYields);
 906     count++;
 907   }
 908
 909   if (fNumberDensitypT)
 910     fNumberDensitypT->Merge(lists[0]);
 911   if (fSumpT)
 912     fSumpT->Merge(lists[1]);
 913   if (fNumberDensityPhi)
 914     fNumberDensityPhi->Merge(lists[2]);
 915   fCorrelationpT->Merge(lists[3]);
 916   fCorrelationEta->Merge(lists[4]);
 917   fCorrelationPhi->Merge(lists[5]);
 918   fCorrelationR->Merge(lists[6]);
 919   fCorrelationLeading2Phi->Merge(lists[7]);
 920   fCorrelationMultiplicity->Merge(lists[8]);
 921   fEventCount->Merge(lists[9]);
 922   fEventCountDifferential->Merge(lists[10]);
 923   fVertexContributors->Merge(lists[11]);
 924   fCentralityDistribution->Merge(lists[12]);
 925   fITSClusterMap->Merge(lists[13]);
 926   if (fTwoTrackDistancePt[0] && lists[14]->GetEntries() > 0)
 927     fTwoTrackDistancePt[0]->Merge(lists[14]);
 928   if (fTwoTrackDistancePt[1] && lists[15]->GetEntries() > 0)
 929     fTwoTrackDistancePt[1]->Merge(lists[15]);
 930   if (fCentralityCorrelation)
 931     fCentralityCorrelation->Merge(lists[16]);
 932   if (fYields && lists[17]->GetEntries() > 0)
 933     fYields->Merge(lists[17]);
 934
 935   for (Int_t i=0; i<kMaxLists; i++)
 936     delete lists[i];
 937
 938   return count+1;
 939 }
 940
 941 void AliUEHistograms::CopyReconstructedData(AliUEHistograms* from)
 942 {
 943   // copies those histograms extracted from ESD to this object
 944
 945   for (Int_t i=0; i<fgkUEHists; i++)
 946     if (GetUEHist(i))
 947       GetUEHist(i)->CopyReconstructedData(from->GetUEHist(i));
 948 }
 949
 950 void AliUEHistograms::DeepCopy(AliUEHistograms* from)
 951 {
 952   // copies the entries of this object's members from the object <from> to this object
 953
 954   for (Int_t i=0; i<fgkUEHists; i++)
 955     if (GetUEHist(i) && from->GetUEHist(i))
 956       GetUEHist(i)->DeepCopy(from->GetUEHist(i));
 957 }
 958
 959 void AliUEHistograms::ExtendTrackingEfficiency(Bool_t verbose)
 960 {
 961   // delegates to AliUEHists
 962
 963   for (Int_t i=0; i<fgkUEHists; i++)
 964     if (GetUEHist(i))
 965       GetUEHist(i)->ExtendTrackingEfficiency(verbose);
 966 }
 967
 968 void AliUEHistograms::Scale(Double_t factor)
 969 {
 970   // scales all contained histograms by the given factor
 971
 972   for (Int_t i=0; i<fgkUEHists; i++)
 973     if (GetUEHist(i))
 974       GetUEHist(i)->Scale(factor);
 975
 976   TList list;
 977   list.Add(fCorrelationpT);
 978   list.Add(fCorrelationEta);
 979   list.Add(fCorrelationPhi);
 980   list.Add(fCorrelationR);
 981   list.Add(fCorrelationLeading2Phi);
 982   list.Add(fCorrelationMultiplicity);
 983   list.Add(fYields);
 984   list.Add(fEventCount);
 985   list.Add(fEventCountDifferential);
 986   list.Add(fVertexContributors);
 987   list.Add(fCentralityDistribution);
 988   list.Add(fCentralityCorrelation);
 989   list.Add(fITSClusterMap);
 990   list.Add(fTwoTrackDistancePt[0]);
 991   list.Add(fTwoTrackDistancePt[1]);
 992
 993   for (Int_t i=0; i<list.GetEntries(); i++)
 994     ((TH1*) list.At(i))->Scale(factor);
 995 }
 996
 997 void AliUEHistograms::Reset()
 998 {
 999   // delegates to AliUEHists
1000
1001   for (Int_t i=0; i<fgkUEHists; i++)
1002     if (GetUEHist(i))
1003       GetUEHist(i)->Reset();
1004 }
1005
1006 Float_t AliUEHistograms::GetInvMassSquared(Float_t pt1, Float_t eta1, Float_t phi1, Float_t pt2, Float_t eta2, Float_t phi2, Float_t m0)
1007 {
1008   // calculate inv mass squared
1009   // same can be achieved, but with more computing time with
1010   /*TLorentzVector photon, p1, p2;
1011   p1.SetPtEtaPhiM(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), 0.510e-3);
1012   p2.SetPtEtaPhiM(particle->Pt(), eta[j], particle->Phi(), 0.510e-3);
1013   photon = p1+p2;
1014   photon.M()*/
1015
1016   Float_t tantheta1 = 1e10;
1017
1018   if (eta1 < -1e-10 || eta1 > 1e-10)
1019     tantheta1 = 2 * TMath::Exp(-eta1) / ( 1 - TMath::Exp(-2*eta1));
1020
1021   Float_t tantheta2 = 1e10;
1022   if (eta2 < -1e-10 || eta2 > 1e-10)
1023     tantheta2 = 2 * TMath::Exp(-eta2) / ( 1 - TMath::Exp(-2*eta2));
1024
1025   Float_t e1squ = m0 * m0 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
1026   Float_t e2squ = m0 * m0 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
1027
1028   Float_t mass2 = 2 * m0 * m0 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( TMath::Cos(phi1 - phi2) + 1.0 / tantheta1 / tantheta2 ) ) );
1029
1030   return mass2;
1031 }