PWGCF/Correlations/DPhi/DiHadronPID/AliAnalysisTaskDiHadronPID.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, 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 // -----------------------------------------------------------------------
  17 //  This analysis task fills histograms with PID information of tracks
  18 //  associated to a high p_T trigger.
  19 // -----------------------------------------------------------------------
  20 //  Author: Misha Veldhoen (misha.veldhoen@cern.ch)
  21
  22 #include <iostream>
  23
  24 // Basic Includes
  25 #include "TH1F.h"
  26 #include "TH2F.h"
  27 #include "TH3F.h"
  28 #include "THn.h"
  29 #include "TFile.h"
  30 #include "TChain.h"
  31 #include "TObject.h"
  32 #include "TObjArray.h"
  33
  34 // Manager/ Handler
  35 #include "AliAnalysisManager.h"
  36 #include "AliInputEventHandler.h"
  37
  38 // Event pool includes.
  39 #include "AliEventPoolManager.h"
  40
  41 // PID includes.
  42 #include "AliPIDResponse.h"
  43
  44 // AOD includes.
  45 #include "AliAODEvent.h"
  46 #include "AliAODTrack.h"
  47 #include "AliAODHandler.h"
  48 #include "AliAODVertex.h"
  49 #include "AliAODInputHandler.h"
  50 #include "AliAODMCParticle.h"
  51 #include "AliAODMCHeader.h"
  52
  53 // Additional includes.
  54 #include "AliTrackDiHadronPID.h"
  55 #include "AliAODTrackCutsDiHadronPID.h"
  56 #include "AliAODEventCutsDiHadronPID.h"
  57 #include "AliHistToolsDiHadronPID.h"
  58 #include "AliFunctionsDiHadronPID.h"
  59
  60 // AnalysisTask headers.
  61 #include "AliAnalysisTaskSE.h"
  62 #include "AliAnalysisTaskDiHadronPID.h"
  63
  64 using namespace std;
  65
  66 ClassImp(AliAnalysisTaskDiHadronPID);
  67
  68 // -----------------------------------------------------------------------
  69 AliAnalysisTaskDiHadronPID::AliAnalysisTaskDiHadronPID():
  70         AliAnalysisTaskSE(),
  71         fPIDResponse(0x0),
  72         fEventCuts(0x0),
  73         fTrackCutsTrigger(0x0),
  74         fTrackCutsAssociated(0x0),
  75         fPoolMgr(0x0),
  76         fTriggerTracks(0x0),
  77         fAssociatedTracks(0x0),
  78         fCurrentAODEvent(0x0),
  79         fOutputList(0x0),
  80         fPtSpectrum(0x0),
  81         fCorrelations(0x0),
  82         fMixedEvents(0x0),
  83         fTOFhistos(0x0),
  84         fTOFPtAxis(0x0),
  85         fTOFTPChistos(0x0),
  86         fTOFTPCPtAxis(0x0),
  87         fNDEtaBins(32),
  88         fNDPhiBins(32),
  89         fMinNEventsForMixing(5),
  90         fPoolTrackDepth(2000),
  91         fPoolSize(1000),
  92         fMixEvents(kTRUE),
  93         fMixTriggers(kFALSE),
  94         fCalculateTOFmismatch(kTRUE),
  95         fT0Fill(0x0),
  96         fLvsEta(0x0),
  97         fLvsEtaProjections(0x0),
  98         fMakeTOFcorrelations(kTRUE),
  99         fMakeTOFTPCcorrelations(kFALSE),
 100         fDebug(0)
 101
 102 {
 103
 104         //
 105         // Default Constructor.
 106         //
 107
 108         if (fDebug > 0) {AliInfo("AliAnalysisTaskDiHadronPID Default Constructor.");}
 109
 110 }
 111
 112 // -----------------------------------------------------------------------
 113 AliAnalysisTaskDiHadronPID::AliAnalysisTaskDiHadronPID(const char* name):
 114         AliAnalysisTaskSE(name),
 115         fPIDResponse(0x0),
 116         fEventCuts(0x0),
 117         fTrackCutsTrigger(0x0),
 118         fTrackCutsAssociated(0x0),
 119         fPoolMgr(0x0),
 120         fTriggerTracks(0x0),
 121         fAssociatedTracks(0x0),
 122         fCurrentAODEvent(0x0),
 123         fOutputList(0x0),
 124         fPtSpectrum(0x0),
 125         fCorrelations(0x0),
 126         fMixedEvents(0x0),
 127         fTOFhistos(0x0),
 128         fTOFPtAxis(0x0),
 129         fTOFTPChistos(0x0),
 130         fTOFTPCPtAxis(0x0),
 131         fNDEtaBins(32),
 132         fNDPhiBins(32),
 133         fMinNEventsForMixing(5),
 134         fPoolTrackDepth(2000),
 135         fPoolSize(1000),
 136         fMixEvents(kTRUE),
 137         fMixTriggers(kFALSE),
 138         fCalculateTOFmismatch(kTRUE),
 139         fT0Fill(0x0),
 140         fLvsEta(0x0),
 141         fLvsEtaProjections(0x0),
 142         fMakeTOFcorrelations(kTRUE),
 143         fMakeTOFTPCcorrelations(kFALSE),
 144         fDebug(0)
 145
 146 {
 147
 148         //
 149         // Named Constructor.
 150         //
 151
 152         if (fDebug > 0) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 153
 154         DefineInput(0,TChain::Class());
 155         DefineOutput(1,TList::Class());
 156
 157 }
 158
 159 // -----------------------------------------------------------------------
 160 AliAnalysisTaskDiHadronPID::~AliAnalysisTaskDiHadronPID() {;
 161
 162         //
 163         // Destructor.
 164         //
 165
 166         if (fDebug > 0) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 167
 168         if (fPoolMgr) {delete fPoolMgr; fPoolMgr = 0x0;}
 169         if (fOutputList) {delete fOutputList; fOutputList = 0x0;}
 170
 171 }
 172
 173 // -----------------------------------------------------------------------
 174 void AliAnalysisTaskDiHadronPID::UserCreateOutputObjects() {
 175
 176         //
 177         // Create Output objects.
 178         //
 179
 180         if (fDebug > 0) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 181
 182         // --- BEGIN: Initialization on the worker nodes ---
 183         AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
 184         AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*> (manager->GetInputEventHandler());
 185
 186         // Getting the pointer to the PID response object.
 187         fPIDResponse = inputHandler->GetPIDResponse();
 188
 189         // Not very neat - only set up for 0-5% analysis.
 190         Int_t nCentralityBins  = 15;
 191         Double_t centralityBins[] = {0., 1., 2., 3., 4., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100.1 };
 192
 193         Int_t nZvtxBins  = 7;
 194         Double_t vertexBins[] = {-7., -5., -3., -1., 1., 3., 5., 7.};
 195
 196         fPoolMgr = new AliEventPoolManager(fPoolSize, fPoolTrackDepth, nCentralityBins, (Double_t*) centralityBins, nZvtxBins, (Double_t*) vertexBins);
 197     // --- END ---
 198
 199         // Create the output list.
 200         fOutputList = new TList();
 201         fOutputList->SetOwner(kTRUE);
 202
 203         // Creating all requested histograms locally.
 204         fEventCuts->CreateHistos();
 205         fOutputList->Add(fEventCuts);
 206
 207         fTrackCutsTrigger->CreateHistos();
 208         fOutputList->Add(fTrackCutsTrigger);
 209
 210         fTrackCutsAssociated->CreateHistos();
 211         fOutputList->Add(fTrackCutsAssociated);
 212
 213         // Get the pT axis for the TOF PID correlations.
 214         Double_t* ptaxis = fTrackCutsAssociated->GetPtAxisPID();
 215         Int_t nptbins = fTrackCutsAssociated->GetNPtBinsPID();
 216         fTOFPtAxis = new TAxis(nptbins, ptaxis);
 217         fTOFPtAxis->SetName("fTOFPtAxis");
 218         fTOFPtAxis->SetTitle("p_{T} GeV/c");
 219         fOutputList->Add(fTOFPtAxis);
 220
 221         // Create Pt spectrum histogram.
 222         fPtSpectrum = new TH1F("fPtSpectrum","p_{T} Spectrum;p_{T} (GeV/c);Count",nptbins,ptaxis);
 223         fOutputList->Add(fPtSpectrum);
 224
 225         // Create unidentified correlations histogram.
 226         fCorrelations = AliHistToolsDiHadronPID::MakeHist3D("fCorrelations","Correlations;#Delta#phi;#Delta#eta;p_{T} (GeV/c)",
 227                 fNDPhiBins,-TMath::Pi()/2.,3.*TMath::Pi()/2.,
 228                 fNDEtaBins,-1.6,1.6,
 229                 nptbins, ptaxis);
 230         fOutputList->Add(fCorrelations);
 231
 232         // Create unidentified mixed events histogram.
 233         fMixedEvents = AliHistToolsDiHadronPID::MakeHist3D("fMixedEvents","Mixed Events;#Delta#phi;#Delta#eta;p_{T} (GeV/c)",
 234                 fNDPhiBins,-TMath::Pi()/2.,3.*TMath::Pi()/2.,
 235                 fNDEtaBins,-1.6,1.6,
 236                 nptbins, ptaxis);
 237         fOutputList->Add(fMixedEvents);
 238
 239         TString speciesname[] = {"Pion","Kaon","Proton"};
 240
 241         // Create TOF correlations histograms (DPhi,DEta,TOF).
 242         if (fMakeTOFcorrelations) {
 243
 244                 fTOFhistos = new TObjArray(3);
 245                 fTOFhistos->SetOwner(kTRUE);
 246                 fTOFhistos->SetName("CorrelationsTOF");
 247
 248                 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
 249
 250                         TObjArray* atmp = new TObjArray(fTOFPtAxis->GetNbins());
 251                         atmp->SetOwner(kTRUE);
 252                         atmp->SetName(speciesname[iSpecies].Data());
 253
 254                         for (Int_t iBinPt = 1; iBinPt < (fTOFPtAxis->GetNbins() + 1); iBinPt++) {
 255
 256                                 Int_t iPtClass = fTrackCutsAssociated->GetPtClass(iBinPt);
 257
 258                                 Int_t NBinsTOF = fTrackCutsAssociated->GetNTOFbins(iPtClass,iSpecies);
 259                                 Double_t TOFmin = fTrackCutsAssociated->GetTOFmin(iPtClass,iSpecies);
 260                                 Double_t TOFmax = fTrackCutsAssociated->GetTOFmax(iPtClass,iSpecies);
 261
 262                                 cout << "ptbin: "<< iBinPt << " class: " << iPtClass << " TOFBins: " << NBinsTOF << " min: " << TOFmin << " max: " << TOFmax << endl;
 263
 264                                 TH3F* htmp = new TH3F(Form("fCorrelationsTOF_%i",iBinPt),
 265                                         Form("%5.3f < p_{T} < %5.3f; #Delta#phi; #Delta#eta; t_{TOF} (ps)", fTOFPtAxis->GetBinLowEdge(iBinPt), fTOFPtAxis->GetBinUpEdge(iBinPt)),
 266                                         fNDPhiBins, -TMath::Pi()/2., 3.*TMath::Pi()/2.,
 267                                         fNDEtaBins, -1.6, 1.6, NBinsTOF, TOFmin, TOFmax);
 268                                 htmp->SetDirectory(0);
 269
 270                                 atmp->Add(htmp);
 271
 272                         }
 273
 274                         fTOFhistos->Add(atmp);
 275
 276                 }
 277
 278                 fOutputList->Add(fTOFhistos);
 279
 280         }
 281
 282         // Create TOF/TPC correlation histograms. (DPhi,DEta,TOF,TPC).
 283         if (fMakeTOFTPCcorrelations) {
 284
 285                 Double_t ptarrayTOFTPC[16] = {2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6,
 286                                                                           2.8, 3.0, 3.2, 3.4, 3.6, 3.8,
 287                                                                           4.2, 4.6, 5.0};
 288                 fTOFTPCPtAxis = new TAxis(15, ptarrayTOFTPC);
 289                 fTOFTPCPtAxis->SetName("fTOFTPCPtAxis");
 290                 fTOFTPCPtAxis->SetTitle("p_{T} GeV/c");
 291                 fOutputList->Add(fTOFTPCPtAxis);
 292
 293                 Int_t NBinsTOFTPC[4] = {32, 32, 60, 40};
 294                 Double_t minTOFTPC[4] = {-TMath::Pi()/2., -1.6, -1., -1.};
 295                 Double_t maxTOFTPC[4] = {3.*TMath::Pi()/2., 1.6, -1., -1.};
 296
 297                 Double_t sTOFest = 110.;
 298                 Double_t sTPCest = 4.5;
 299
 300                 fTOFTPChistos = new TObjArray(3);
 301                 fTOFTPChistos->SetOwner(kTRUE);
 302                 fTOFTPChistos->SetName("TOFTPChistos");
 303
 304                 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
 305
 306                         TObjArray* atmp = new TObjArray(fTOFTPCPtAxis->GetNbins());
 307                         atmp->SetOwner(kTRUE);
 308                         atmp->SetName(speciesname[iSpecies].Data());
 309
 310                         for (Int_t iBinPt = 1; iBinPt < (fTOFTPCPtAxis->GetNbins() + 1); iBinPt++) {
 311
 312                                 // Determine PID ranges.
 313
 314                                 // Set range +/- 5 sigma of main peak. (+ 10 sigma for TOF max, for mismatches.)
 315                                 Double_t TOFmin = -5. * sTOFest;
 316                                 Double_t TOFmax = 10. * sTOFest;
 317                                 Double_t TPCmin = -4. * sTPCest;
 318                                 Double_t TPCmax = 4. * sTPCest;
 319
 320                                 Double_t TOFexp = AliFunctionsDiHadronPID::TOFExpTime(fTOFTPCPtAxis->GetBinLowEdge(iBinPt), 0.4,  AliFunctionsDiHadronPID::M(iSpecies));
 321                                 Double_t TPCexp = AliFunctionsDiHadronPID::TPCExpdEdX(fTOFTPCPtAxis->GetBinLowEdge(iBinPt), 0.4,  AliFunctionsDiHadronPID::M(iSpecies));
 322
 323                                 for (Int_t jSpecies = 0; jSpecies < 3; jSpecies++) {
 324
 325                                         if (iSpecies == jSpecies) {continue;}
 326
 327                                         Double_t TOFexpOther = AliFunctionsDiHadronPID::TOFExpTime(fTOFTPCPtAxis->GetBinLowEdge(iBinPt), 0.4,  AliFunctionsDiHadronPID::M(jSpecies));
 328                                         Double_t TPCexpOther = AliFunctionsDiHadronPID::TPCExpdEdX(fTOFTPCPtAxis->GetBinLowEdge(iBinPt), 0.4,  AliFunctionsDiHadronPID::M(jSpecies));
 329
 330                                         // If any peak is within +/- 7 sigma, then also add this peak.
 331                                         if ( (TMath::Abs(TOFexp - TOFexpOther) < 7. * sTOFest) ||
 332                                                  (TMath::Abs(TPCexp - TPCexpOther) < 7. * sTPCest) ) {
 333
 334                                                 TOFmin = TMath::Min(TOFmin, (TOFexpOther - TOFexp - 2. * sTOFest) );
 335                                                 TOFmax = TMath::Max(TOFmax, (TOFexpOther - TOFexp + 10. * sTOFest) );
 336                                                 TPCmin = TMath::Min(TPCmin, (TPCexpOther - TPCexp - 2. * sTPCest) );
 337                                                 TPCmax = TMath::Max(TPCmax, (TPCexpOther - TPCexp + 2. * sTPCest) );
 338
 339                                         }
 340
 341                                 }
 342
 343                                 minTOFTPC[2] = TOFmin;
 344                                 maxTOFTPC[2] = TOFmax;
 345                                 minTOFTPC[3] = TPCmin;
 346                                 maxTOFTPC[3] = TPCmax;
 347
 348                                 THnF* htmp = new THnF(Form("fCorrelationsTOF_%i",iBinPt),
 349                                         Form("%5.3f < p_{T} < %5.3f", fTOFTPCPtAxis->GetBinLowEdge(iBinPt), fTOFTPCPtAxis->GetBinUpEdge(iBinPt)),
 350                                         4, NBinsTOFTPC, minTOFTPC, maxTOFTPC);
 351
 352                                 (htmp->GetAxis(0))->SetTitle("#Delta#phi");
 353                                 (htmp->GetAxis(1))->SetTitle("#Delta#eta");
 354                                 (htmp->GetAxis(2))->SetTitle("t_{TOF} (ps)");
 355                                 (htmp->GetAxis(3))->SetTitle("dE/dx (a.u.)");
 356
 357                                 atmp->Add(htmp);
 358
 359                         }
 360
 361                         fTOFTPChistos->Add(atmp);
 362
 363                 }
 364
 365                 fOutputList->Add(fTOFTPChistos);
 366
 367         }
 368
 369         // Load external TOF histograms if flag is set.
 370         if (fCalculateTOFmismatch) {LoadExtMismatchHistos();}
 371
 372         PostData(1,fOutputList);
 373
 374 }
 375
 376 // -----------------------------------------------------------------------
 377 void AliAnalysisTaskDiHadronPID::LocalInit() {
 378
 379         //
 380         // Initialize on the this computer.
 381         //
 382
 383         if (fDebug > 0) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 384
 385 }
 386
 387 // -----------------------------------------------------------------------
 388 void AliAnalysisTaskDiHadronPID::UserExec(Option_t*) {
 389
 390         //
 391         // Main Loop.
 392         //
 393
 394         if (fDebug > 0) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 395
 396         // Input Current Event.
 397         fCurrentAODEvent = dynamic_cast<AliAODEvent*>(InputEvent());
 398         if (!fCurrentAODEvent) AliFatal("No Event Found!");
 399
 400         if (!fEventCuts->IsSelected(fCurrentAODEvent)) {return;}
 401
 402         // Fill the global tracks array. - NOT NEEDED I THINK, since we're not using
 403         // bit 1<<7 for the associated tracks!
 404
 405         // Let the track cut objects know that a new event will start.
 406         fTrackCutsTrigger->StartNewEvent();
 407         fTrackCutsAssociated->StartNewEvent();
 408
 409         // Create arrays for trigger/associated particles.
 410         fTriggerTracks = new TObjArray(350);
 411         fTriggerTracks->SetOwner(kTRUE);
 412
 413         fAssociatedTracks = new TObjArray(3500);
 414         fAssociatedTracks->SetOwner(kTRUE);
 415
 416         for (Int_t iTrack = 0; iTrack < fCurrentAODEvent->GetNumberOfTracks(); iTrack++) {
 417
 418                 AliAODTrack* track = (AliAODTrack*)fCurrentAODEvent->GetTrack(iTrack);
 419                 AliTrackDiHadronPID* pidtrack = new AliTrackDiHadronPID(track,0x0,0x0,fPIDResponse);
 420                 pidtrack->ForgetAboutPointers();
 421                 pidtrack->SetDebugLevel(fDebug);
 422
 423                 Double_t rndhittime = -1.e21;
 424                 if (fCalculateTOFmismatch) rndhittime = GenerateRandomHit(pidtrack->Eta());
 425
 426                 // Fill p_T spectrum.
 427                 fPtSpectrum->Fill(pidtrack->Pt());
 428
 429                 // Fill the trigger/associated tracks array.
 430                 if (fTrackCutsTrigger->IsSelectedData(pidtrack,rndhittime)) {fTriggerTracks->AddLast(pidtrack);}
 431                 else if (fTrackCutsAssociated->IsSelectedData(pidtrack,rndhittime)) {fAssociatedTracks->AddLast(pidtrack);}
 432                 else {delete pidtrack;}
 433
 434         }
 435
 436         // Fill Correlation histograms.
 437         for (Int_t iTrigger = 0; iTrigger < fTriggerTracks->GetEntriesFast(); iTrigger++) {
 438                 AliTrackDiHadronPID* triggertrack = (AliTrackDiHadronPID*)fTriggerTracks->At(iTrigger);
 439
 440                 for (Int_t iAssociated = 0; iAssociated < fAssociatedTracks->GetEntriesFast(); iAssociated++) {
 441                         AliTrackDiHadronPID* associatedtrack = (AliTrackDiHadronPID*)fAssociatedTracks->At(iAssociated);
 442
 443                         Double_t DPhi = triggertrack->Phi() - associatedtrack->Phi();
 444                         if (DPhi < -TMath::Pi()/2.) {DPhi += 2.*TMath::Pi();}
 445                         if (DPhi > 3.*TMath::Pi()/2.) {DPhi -= 2.*TMath::Pi();}
 446
 447                         Double_t DEta = triggertrack->Eta() - associatedtrack->Eta();
 448                         fCorrelations->Fill(DPhi,DEta,associatedtrack->Pt());
 449
 450                         for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
 451
 452                                 Double_t apt = associatedtrack->Pt();
 453
 454                                 // Fill TOF correlations.
 455                                 if (fMakeTOFcorrelations) {
 456
 457                                         TObjArray* atmp = (TObjArray*)fTOFhistos->At(iSpecies);
 458                                         Int_t ptbin = fTOFPtAxis->FindBin(apt);
 459
 460                                         // Only fill if histogram exists in fTOFhistos.
 461                                         if ( !(ptbin < 1) && !(ptbin > fTOFPtAxis->GetNbins()) ) {
 462
 463                                                 TH3F* htmp = (TH3F*)atmp->At(ptbin - 1);
 464                                                 htmp->Fill(DPhi, DEta, associatedtrack->GetTOFsignalMinusExpected(iSpecies));
 465
 466                                         }
 467                                 }
 468
 469                                 // Fill TOF/ TPC Correlations.
 470                                 if (fMakeTOFTPCcorrelations) {
 471
 472                                         TObjArray* atmp = (TObjArray*)fTOFTPChistos->At(iSpecies);
 473                                         Int_t ptbin = fTOFTPCPtAxis->FindBin(apt);
 474
 475                                         // Only fill if histogram exists in fTOFhistos.
 476                                         if ( !(ptbin < 1) && !(ptbin > fTOFTPCPtAxis->GetNbins()) ) {
 477
 478                                                 THnF* htmp = (THnF*)atmp->At(ptbin - 1);
 479                                                 Double_t TOFTPCfill[4] = {DPhi, DEta,
 480                                                         associatedtrack->GetTOFsignalMinusExpected(iSpecies), associatedtrack->GetTPCsignalMinusExpected(iSpecies)};
 481
 482                                                 htmp->Fill(TOFTPCfill);
 483
 484                                         }
 485                                 }
 486                         }
 487                 }
 488         }
 489
 490         cout<<"Triggers: "<<fTriggerTracks->GetEntriesFast()<<" Associateds: "<<fAssociatedTracks->GetEntriesFast()<<endl;
 491
 492         // Determine centrality of current event.
 493         TString centralityestimator = fEventCuts->GetCentralityEstimator();
 494         AliCentrality* currentcentrality = fCurrentAODEvent->GetCentrality();
 495         Float_t percentile = currentcentrality->GetCentralityPercentile(centralityestimator.Data());
 496
 497         // Determine vtxz of current event.
 498         AliAODVertex* currentprimaryvertex = fCurrentAODEvent->GetPrimaryVertex();
 499         Double_t vtxz = currentprimaryvertex->GetZ();
 500
 501         AliEventPool* poolin = fPoolMgr->GetEventPool(percentile, vtxz);
 502         if (!poolin) {AliFatal(Form("No pool found for centrality = %f, vtxz = %f", percentile, vtxz));}
 503         // TObjArray* fGlobalTracksArray;
 504
 505         // Give a print out of the pool manager's contents.
 506         if (fDebug > 0) PrintPoolManagerContents();
 507
 508         // Mix events if there are enough events in the pool.
 509         if (poolin->GetCurrentNEvents() >= fMinNEventsForMixing) {
 510                 {cout << "Mixing Events." << endl;}
 511
 512                 // Loop over all events in the event pool.
 513                 for (Int_t iMixEvent = 0; iMixEvent < poolin->GetCurrentNEvents(); iMixEvent++) {
 514                 TObjArray* mixtracks = poolin->GetEvent(iMixEvent);
 515
 516                 // Mix either the triggers or the associateds.
 517                 if (fMixTriggers) {
 518
 519                                 // Loop over all associateds in this event.
 520                                 for (Int_t iAssociated = 0; iAssociated < fAssociatedTracks->GetEntriesFast(); iAssociated++) {
 521                                         AliTrackDiHadronPID* associatedtrack = (AliTrackDiHadronPID*)fAssociatedTracks->At(iAssociated);
 522
 523                                 // Loop over all mixed tracks.
 524                                 for (Int_t iMixTrack = 0; iMixTrack < mixtracks->GetEntriesFast(); iMixTrack++) {
 525                                         AliTrackDiHadronPID* mixtrack = (AliTrackDiHadronPID*)mixtracks->At(iMixTrack);
 526
 527                                                 Double_t DPhi = mixtrack->Phi() - associatedtrack->Phi();
 528                                                 if (DPhi < -TMath::Pi()/2.) {DPhi += 2.*TMath::Pi();}
 529                                                 if (DPhi > 3.*TMath::Pi()/2.) {DPhi -= 2.*TMath::Pi();}
 530
 531                                                 Double_t DEta = mixtrack->Eta() - associatedtrack->Eta();
 532                                                 fMixedEvents->Fill(DPhi,DEta,associatedtrack->Pt());
 533
 534                                 }
 535                                 }
 536
 537                         } else {
 538
 539                                 // Loop over all triggers in this event.
 540                                 for (Int_t iTrigger = 0; iTrigger < fTriggerTracks->GetEntriesFast(); iTrigger++) {
 541                                         AliTrackDiHadronPID* triggertrack = (AliTrackDiHadronPID*)fTriggerTracks->At(iTrigger);
 542
 543                                 // Loop over all mixed tracks.
 544                                 for (Int_t iMixTrack = 0; iMixTrack < mixtracks->GetEntriesFast(); iMixTrack++) {
 545                                         AliTrackDiHadronPID* mixtrack = (AliTrackDiHadronPID*)mixtracks->At(iMixTrack);
 546
 547                                                 Double_t DPhi = triggertrack->Phi() - mixtrack->Phi();
 548                                                 if (DPhi < -TMath::Pi()/2.) {DPhi += 2.*TMath::Pi();}
 549                                                 if (DPhi > 3.*TMath::Pi()/2.) {DPhi -= 2.*TMath::Pi();}
 550
 551                                                 Double_t DEta = triggertrack->Eta() - mixtrack->Eta();
 552                                                 fMixedEvents->Fill(DPhi,DEta,mixtrack->Pt());
 553
 554                                 }
 555                                 }
 556
 557                         } // End if
 558                 }
 559         }
 560
 561         // Update the event pool.
 562         AliEventPool* poolout = fPoolMgr->GetEventPool(percentile, vtxz); // Get the buffer associated with the current centrality and z-vtx
 563         if (!poolout) AliFatal(Form("No pool found for centrality = %f, vtx_z = %f", percentile, vtxz));
 564
 565         // Q: is it a problem that the fAssociatedTracks array can be bigger than the number of tracks inside?
 566         if (fMixTriggers) {
 567                 poolout->UpdatePool(fTriggerTracks);
 568                 fAssociatedTracks->Delete();
 569                 delete fAssociatedTracks;
 570         }
 571         else {
 572                 poolout->UpdatePool(fAssociatedTracks);
 573                 fTriggerTracks->Delete();
 574                 delete fTriggerTracks;
 575         }
 576
 577         fTriggerTracks = 0x0;
 578         fAssociatedTracks = 0x0;
 579
 580         // Tell the track cut object that the event is done.
 581         fTrackCutsTrigger->EventIsDone(0);
 582         fTrackCutsAssociated->EventIsDone(0);
 583
 584         PostData(1,fOutputList);
 585
 586 }
 587
 588 // -----------------------------------------------------------------------
 589 Bool_t AliAnalysisTaskDiHadronPID::LoadExtMismatchHistos() {
 590
 591         //
 592         // Attempting to load a root file containing information needed
 593         // to generate random TOF hits.
 594         //
 595
 596         if (fDebug > 0) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 597
 598         // Opening external TOF file.
 599         if (fDebug > 0) cout<<"Trying to open TOFmismatchHistos.root ..."<<endl;
 600         TFile* fin = 0x0;
 601         fin = TFile::Open("alien:///alice/cern.ch/user/m/mveldhoe/rootfiles/TOFmismatchHistos.root");
 602         if (!fin) {
 603                 AliWarning("Couln't open TOFmismatchHistos, will not calculate mismatches...");
 604                 fCalculateTOFmismatch = kFALSE;
 605                 return kFALSE;
 606         }
 607
 608         // Check if the required histograms are present.
 609         TH1F* tmp1 = (TH1F*)fin->Get("hNewT0Fill");
 610         if (!tmp1) {
 611                 AliWarning("Couln't find hNewT0Fill, will not calculate mismatches...");
 612                 fCalculateTOFmismatch = kFALSE;
 613                 return kFALSE;
 614         }
 615         TH2F* tmp2 = (TH2F*)fin->Get("hLvsEta");
 616         if (!tmp2) {
 617                 AliWarning("Couln't find hLvsEta, will not calculate mismatches...");
 618                 fCalculateTOFmismatch = kFALSE;
 619                 return kFALSE;
 620         }
 621
 622         // Make a deep copy of the files in the histogram.
 623         fT0Fill = (TH1F*)tmp1->Clone("fT0Fill");
 624         fLvsEta = (TH2F*)tmp2->Clone("fLvsEta");
 625
 626         // Close the external file.
 627         AliInfo("Closing external file.");
 628         fin->Close();
 629
 630         // Creating a TObjArray for LvsEta projections.
 631         const Int_t nbinseta = fLvsEta->GetNbinsX();
 632         fLvsEtaProjections = new TObjArray(nbinseta);
 633         fLvsEtaProjections->SetOwner(kTRUE);
 634
 635         // Making the projections needed (excluding underflow/ overflow).
 636         for (Int_t iEtaBin = 1; iEtaBin < (nbinseta + 1); iEtaBin++) {
 637                 TH1F* tmp = (TH1F*)fLvsEta->ProjectionY(Form("LvsEtaProjection_%i",iEtaBin),iEtaBin,iEtaBin);
 638                 tmp->SetDirectory(0);
 639                 fLvsEtaProjections->AddAt(tmp,iEtaBin);
 640         }
 641
 642         return kTRUE;
 643
 644 }
 645
 646 // -----------------------------------------------------------------------
 647 Double_t AliAnalysisTaskDiHadronPID::GenerateRandomHit(Double_t eta) {
 648
 649         //
 650         // Returns a random TOF time.
 651         //
 652
 653         if (fDebug > 0) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 654
 655         // Default (error) value:
 656         Double_t rndhittime = -1.e21;
 657
 658         // TOF mismatch flag is not turned on.
 659         if (!fCalculateTOFmismatch) {
 660                 AliFatal("Called GenerateRandomHit() method, but flag fCalculateTOFmismatch not set.");
 661                 return rndhittime;
 662         }
 663
 664         // TOF doesn't extend much further than 0.8.
 665         if (TMath::Abs(eta) > 0.8) {
 666                 if (fDebug) {AliInfo("Tried to get a random hit for a track with eta > 0.8.");}
 667                 return rndhittime;
 668         }
 669
 670         // Finding the bin of the eta.
 671         TAxis* etaAxis = fLvsEta->GetXaxis();
 672         Int_t etaBin = etaAxis->FindBin(eta);
 673         //cout<<"Eta: "<<eta<<" bin: "<<etaBin<<endl;
 674         const TH1F* lengthDistribution = (const TH1F*)fLvsEtaProjections->At(etaBin);
 675
 676         if (!lengthDistribution) {
 677                 AliFatal("length Distribution not found.");
 678                 return rndhittime;
 679         }
 680
 681         Double_t currentRndLength = lengthDistribution->GetRandom(); // in cm.
 682
 683         // Similar to Roberto's code.
 684         Double_t currentRndTime = currentRndLength / (TMath::C() * 1.e2 / 1.e12);
 685         Double_t t0fill = -1.26416e+04;
 686         rndhittime = fT0Fill->GetRandom() - t0fill + currentRndTime;
 687
 688         return rndhittime;
 689
 690 }
 691
 692 // -----------------------------------------------------------------------
 693 void AliAnalysisTaskDiHadronPID::PrintPoolManagerContents() {
 694
 695         //
 696         // Prints out the current contents of the event pool manager.
 697         //
 698
 699         // Determine the number of pools in the pool manager.
 700         AliEventPool* poolin = fPoolMgr->GetEventPool(0,0);
 701         Int_t NPoolsCentrality = 0;
 702         while (poolin) {
 703                 NPoolsCentrality++;
 704                 poolin = fPoolMgr->GetEventPool(NPoolsCentrality,0);
 705         }
 706
 707         poolin = fPoolMgr->GetEventPool(0,0);
 708         Int_t NPoolsVtxZ = 0;
 709         while (poolin) {
 710                 NPoolsVtxZ++;
 711                 poolin = fPoolMgr->GetEventPool(0,NPoolsVtxZ);
 712         }
 713
 714         // Loop over all Pools in the matrix of the pool manager.
 715         cout<<" Pool manager contents: (Nevt,NTrack)"<<endl;
 716         for (Int_t iCentrality = 0; iCentrality < NPoolsCentrality; iCentrality++) {
 717                 cout<<Form("Centrality Bin: %2i --> ", iCentrality);
 718
 719                 for (Int_t iVtxZ = 0; iVtxZ < NPoolsVtxZ; iVtxZ++) {
 720
 721                         poolin = fPoolMgr->GetEventPool(iCentrality, iVtxZ);
 722
 723                         cout<<Form("(%2i,%4i) ",poolin->GetCurrentNEvents(), poolin->NTracksInPool());
 724
 725                 }
 726
 727                 cout<<endl;
 728         }
 729
 730 }
 731
 732 // -----------------------------------------------------------------------
 733 void AliAnalysisTaskDiHadronPID::Terminate(Option_t*) {;
 734
 735         //
 736         // Called when task is done.
 737         //
 738
 739         if (fDebug > 0) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
 740
 741         delete fT0Fill;
 742         fT0Fill = 0x0;
 743         delete fLvsEta;
 744         fLvsEta = 0x0;
 745         delete fLvsEtaProjections;
 746         fLvsEtaProjections = 0x0;
 747
 748 }