PWGCF/Correlations/DPhi/AliAnalysisTaskDiHadronPID.cxx

   1 // ----------------------------------------------------------------------------
   2 // This class makes di-hadron correlations, with TOF and TPC signals for
   3 // the associated particles.
   4 //
   5 // ----------------------------------------------------------------------------
   6 // Author: Misha Veldhoen (misha.veldhoen@cern.ch)
   7 // Last edit: May 2nd 2012. (v 8.00)
   8 // ----------------------------------------------------------------------------
   9
  10 #include <iostream>
  11 #include "TChain.h"
  12 #include "TTree.h"
  13 #include "TH1F.h"
  14 #include "TH2F.h"
  15 #include "TH3F.h"
  16 #include "TCanvas.h"
  17 #include "TFile.h"
  18 #include "TClonesArray.h"
  19
  20 #include "AliAODTrack.h"
  21 #include "AliAODEvent.h"
  22 #include "AliAODInputHandler.h"
  23 #include "AliAODVertex.h"
  24 //#include "AliAODPid.h"
  25
  26 #include "AliAnalysisManager.h"
  27 #include "AliInputEventHandler.h"
  28 #include "AliPIDResponse.h"
  29 #include "AliTPCPIDResponse.h"
  30 //#include "AliTOFPIDResponse.h"
  31
  32 // Includes for a MC run.
  33 #include "AliAODMCParticle.h"
  34
  35 using namespace std;
  36
  37 #include "AliAnalysisTaskDiHadronPID.h"
  38
  39 ClassImp(AliAnalysisTaskDiHadronPID);
  40
  41 //_____________________________________________________________________________
  42 AliAnalysisTaskDiHadronPID::AliAnalysisTaskDiHadronPID():
  43         AliAnalysisTaskSE(),
  44         fPIDResponse(0x0),
  45         fAODEvent(0x0),
  46         fAODHeader(0x0),
  47         fAODVertex(0x0),
  48         fAODTrack(0x0),
  49         fGlobalTracks(0x0),
  50         fMCTracks(0x0),
  51         fCentrality(0x0),
  52         fVertexZ(0x0),
  53     fDCA(0x0),
  54     fDCAZoomed(0x0),
  55     fDCAZoomedTwice(0x0),
  56     fDCACut(0x0),
  57     fDCAZoomedCut(0x0),
  58     fDCAZoomedTwiceCut(0x0),
  59     fITSHits(0x0),
  60     fTrackCutsCount(0x0),
  61     fTrackCutsPt(0x0),
  62     fTrackCutsEta(0x0),
  63     fTrackCutsPhi(0x0),
  64     fEtaSpectrumTrig(0x0),
  65     fEtaSpectrumAssoc(0x0),
  66     fPhiSpectrumAssoc(0x0),
  67         fTPCnSigmaProton(0x0),
  68         fTOFnSigmaProton(0x0),
  69         fTPCnSigmaPion(0x0),
  70         fTOFnSigmaPion(0x0),
  71         fTPCnSigmaKaon(0x0),
  72         fTOFnSigmaKaon(0x0),
  73         fTPCSignal(0x0),
  74         fTOFSignal(0x0),
  75         fDiHadron(0x0),
  76     fMixedEvents(0x0),
  77         fHistoList(0x0),
  78     fCalculateMixedEvents(kFALSE),
  79     fBeamType("PbPb"),
  80     fMC(kFALSE),
  81     fMaxEta(0.8),
  82     fMaxPlotEta(0.8),
  83     fMaxRap(0.5),
  84     fMaxPt(10.),
  85     fNEtaBins(32),
  86     fNPhiBins(36),
  87     fVertexZMixedEvents(2.),
  88     fCentralityCutMax(0.),
  89     fCentralityCutMin(10.),
  90     fZoomed(kFALSE),
  91     fDoITSCut(kFALSE),
  92     fDoDCACut(kFALSE),
  93     fDemandNoMismatch(kFALSE),
  94     fVerbose(0),
  95     fPrintBufferSize(kFALSE),
  96         fTrigBufferIndex(0),
  97         fTrigBufferSize(0),
  98         fTrigBufferMaxSize(100)
  99
 100 {
 101
 102         //
 103         // Default Constructor.
 104         //
 105
 106         // Trigger buffer.
 107         for(Int_t ii=0; ii<25000; ii++) {
 108                 for(Int_t jj=0; jj<4; jj++) {
 109                         fTrigBuffer[ii][jj]=0;
 110                 }
 111         }
 112
 113     // The identified di-hadron correlations.
 114         for (Int_t ii = 0; ii < 3; ii++) {
 115                 //fMixedEventsTPCTOFCut[ii]=0x0;
 116                 for (Int_t jj = 0; jj < 10; jj++) {
 117                         fDiHadronTPCTOF[ii][jj]=0x0;
 118                         fMixedEventsTPCTOF[ii][jj]=0x0;
 119             fInclusiveTPCTOF[ii][jj]=0x0;
 120             fInclusiveTPCTOFRap[ii][jj]=0x0;
 121                 }
 122         }
 123
 124     // Track cut labels.
 125     for (Int_t ii=0; ii<8; ii++) {
 126         fTrackCutLabelNumbers[ii]=ii+1;
 127     }
 128
 129         // Monte Carlo Efficiency Histograms.
 130         for (Int_t iSpecies=0; iSpecies<6; iSpecies++) {
 131         fPtEtaDistrDataPrim[iSpecies]=0x0;
 132         fPtEtaDistrDataSec[iSpecies]=0x0;
 133         fPtRapDistrDataPrimRapCut[iSpecies]=0x0;
 134         fPtRapDistrDataSecRapCut[iSpecies]=0x0;
 135
 136                 fPtEtaDistrMCPrim[iSpecies]=0x0;
 137                 fPtEtaDistrMCSec[iSpecies]=0x0;
 138                 fPtRapDistrMCPrimRapCut[iSpecies]=0x0;
 139                 fPtRapDistrMCSecRapCut[iSpecies]=0x0;
 140
 141                 fDiHadronMC[iSpecies]=0x0;
 142
 143         }
 144
 145 }
 146
 147 //_____________________________________________________________________________
 148 AliAnalysisTaskDiHadronPID::AliAnalysisTaskDiHadronPID(const char *name):
 149         AliAnalysisTaskSE(name),
 150     fPIDResponse(0x0),
 151     fAODEvent(0x0),
 152     fAODHeader(0x0),
 153     fAODVertex(0x0),
 154     fAODTrack(0x0),
 155     fGlobalTracks(0x0),
 156         fMCTracks(0x0),
 157     fCentrality(0x0),
 158     fVertexZ(0x0),
 159     fDCA(0x0),
 160     fDCAZoomed(0x0),
 161     fDCAZoomedTwice(0x0),
 162     fDCACut(0x0),
 163     fDCAZoomedCut(0x0),
 164     fDCAZoomedTwiceCut(0x0),
 165     fITSHits(0x0),
 166     fTrackCutsCount(0x0),
 167     fTrackCutsPt(0x0),
 168     fTrackCutsEta(0x0),
 169     fTrackCutsPhi(0x0),
 170     fEtaSpectrumTrig(0x0),
 171     fEtaSpectrumAssoc(0x0),
 172     fPhiSpectrumAssoc(0x0),
 173     fTPCnSigmaProton(0x0),
 174     fTOFnSigmaProton(0x0),
 175     fTPCnSigmaPion(0x0),
 176     fTOFnSigmaPion(0x0),
 177     fTPCnSigmaKaon(0x0),
 178     fTOFnSigmaKaon(0x0),
 179     fTPCSignal(0x0),
 180     fTOFSignal(0x0),
 181         fDiHadron(0x0),
 182     fMixedEvents(0x0),
 183     fHistoList(0x0),
 184     fCalculateMixedEvents(kFALSE),
 185     fBeamType("PbPb"),
 186     fMC(kFALSE),
 187     fMaxEta(0.8),
 188     fMaxPlotEta(0.8),
 189     fMaxRap(0.5),
 190     fMaxPt(10.),
 191     fNEtaBins(32),
 192     fNPhiBins(36),
 193     fVertexZMixedEvents(2.),
 194     fCentralityCutMax(0.),
 195     fCentralityCutMin(10.),
 196     fZoomed(kFALSE),
 197     fDoITSCut(kFALSE),
 198     fDoDCACut(kFALSE),
 199     fDemandNoMismatch(kFALSE),
 200     fVerbose(0),
 201     fPrintBufferSize(kFALSE),
 202     fTrigBufferIndex(0),
 203     fTrigBufferSize(0),
 204         fTrigBufferMaxSize(100)
 205
 206 {
 207
 208         //
 209         // Named Constructor.
 210         //
 211
 212         DefineInput(0, TChain::Class());
 213         DefineOutput(1, TList::Class());
 214
 215         // Trigger buffer.
 216         for(Int_t ii=0; ii<25000; ii++) {
 217                 for(Int_t jj=0; jj<4; jj++) {
 218                         fTrigBuffer[ii][jj]=0;
 219                 }
 220         }
 221
 222     // The identified di-hadron correlations.
 223         for (Int_t ii = 0; ii < 3; ii++) {
 224                 //fMixedEventsTPCTOFCut[ii]=0x0;
 225                 for (Int_t jj = 0; jj < 10; jj++) {
 226                         fDiHadronTPCTOF[ii][jj]=0x0;
 227                         fMixedEventsTPCTOF[ii][jj]=0x0;
 228             fInclusiveTPCTOF[ii][jj]=0x0;
 229             fInclusiveTPCTOFRap[ii][jj]=0x0;
 230                 }
 231         }
 232
 233     // Track cut labels.
 234     for (Int_t ii=0; ii<8; ii++) {
 235         fTrackCutLabelNumbers[ii]=ii+1;
 236     }
 237
 238         // Monte Carlo Efficiency Histograms.
 239         for (Int_t iSpecies=0; iSpecies<6; iSpecies++) {
 240         fPtEtaDistrDataPrim[iSpecies]=0x0;
 241         fPtEtaDistrDataSec[iSpecies]=0x0;
 242         fPtRapDistrDataPrimRapCut[iSpecies]=0x0;
 243         fPtRapDistrDataSecRapCut[iSpecies]=0x0;
 244
 245                 fPtEtaDistrMCPrim[iSpecies]=0x0;
 246                 fPtEtaDistrMCSec[iSpecies]=0x0;
 247                 fPtRapDistrMCPrimRapCut[iSpecies]=0x0;
 248                 fPtRapDistrMCSecRapCut[iSpecies]=0x0;
 249
 250                 fDiHadronMC[iSpecies]=0x0;
 251
 252         }
 253
 254 }
 255
 256 //_____________________________________________________________________________
 257 AliAnalysisTaskDiHadronPID::~AliAnalysisTaskDiHadronPID() {
 258
 259         //
 260         // Destructor.
 261         //
 262
 263     if(fGlobalTracks) {
 264         delete fGlobalTracks;
 265         fGlobalTracks=0x0;
 266     }
 267         /*
 268         if (fMCTracks) {
 269                 delete fMCTracks;
 270                 fMCTracks=0x0;
 271         }
 272 */
 273 }
 274
 275 //_____________________________________________________________________________
 276 void AliAnalysisTaskDiHadronPID::UserCreateOutputObjects()
 277
 278 {
 279         //
 280         // Output objects.
 281         //
 282
 283     // Print the settings of the analysis task.
 284     cout<<endl;
 285     cout<<"-----------------------------------------------"<<endl;
 286     cout<<" AliAnalysisTaskDiHadronPID Settings:"<<endl;
 287     cout<<"-----------------------------------------------"<<endl;
 288     cout<<"Verbose Level: "<<fVerbose<<endl;
 289     cout<<"Mixed Events Calculated: "<<fCalculateMixedEvents<<endl;
 290     cout<<"Beam Type: "<<fBeamType<<endl;
 291     cout<<"Run over MC: "<<fMC<<endl;
 292     cout<<Form("Max eta: %3.1f",fMaxEta)<<endl;
 293     cout<<Form("Max eta plotted: %3.1f",fMaxPlotEta)<<endl;
 294     cout<<Form("Max rapidity in spectra: %3.1f",fMaxRap)<<endl;
 295     cout<<Form("Max p_T for the triggers: %3.1f GeV/c.",fMaxPt)<<endl;
 296     cout<<"Nbins in eta and delta eta: "<<fNEtaBins<<endl;
 297     cout<<"Nbins in phi and delta phi: "<<fNPhiBins<<endl;
 298     cout<<Form("Events mixed if vertices differ %3.1f cm.",fVertexZMixedEvents)<<endl;
 299     cout<<Form("Centrality between %3.1f and %3.1f percent.",fCentralityCutMax,fCentralityCutMin)<<endl;
 300     cout<<"Tracks are cut if less than 2 SPD hits: "<<fDoITSCut<<endl;
 301     cout<<"Tracks cut if DCA is too big: "<<fDoDCACut<<endl;
 302     cout<<"Tracks cut if there is a TPC-TOF mismatch: "<<fDemandNoMismatch<<endl;
 303     cout<<"Maximum number of triggers stored: "<<fTrigBufferMaxSize<<endl;
 304     cout<<"-----------------------------------------------"<<endl;
 305     cout<<endl;
 306
 307         // Obtain a pointer to the analysis manager.
 308         AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
 309     if (!manager) {
 310         if (fVerbose>0) cout<<"AliAnalysisTaskDiHadronPID::UserCreateOutputObjects -> ERROR: Analysis manager not found."<<endl;
 311         return;
 312     }
 313     if (fVerbose>1) {
 314         cout<<"AliAnalysisTaskDiHadronPID::UserCreateOutputObjects -> Analysis manager found."<<endl;
 315     }
 316
 317     // Obtain a pointer to the input handler.
 318     AliInputEventHandler* inputHandler = dynamic_cast<AliInputEventHandler*> (manager->GetInputEventHandler());
 319     if (!inputHandler) {
 320         if (fVerbose>0) cout<<"AliAnalysisTaskDiHadronPID::UserCreateOutputObjects -> ERROR: Input handler not found."<<endl;
 321         return;
 322     }
 323     if (fVerbose>1) {
 324         cout<<"AliAnalysisTaskDiHadronPID::UserCreateOutputObjects -> Input handler found."<<endl;
 325     }
 326
 327         // Obtain a pointer to the PID response object.
 328         fPIDResponse = inputHandler->GetPIDResponse();
 329     if (!fPIDResponse) {
 330         if (fVerbose>0) cout<<"AliAnalysisTaskDiHadronPID::UserCreateOutputObjects -> ERROR: PID response object not found."<<endl;
 331         return;
 332     }
 333     if (fVerbose>1) {
 334         cout<<"AliAnalysisTaskDiHadronPID::UserCreateOutputObjects -> PID response object found."<<endl;
 335     }
 336
 337         // Create the output of the task.
 338         fHistoList = new TList();
 339         fHistoList->SetOwner(kTRUE);
 340
 341         // Ranges in dPhi, dEta, and the number of bins for di-hadron correlations.
 342         Int_t binsHisto[4] = {fNPhiBins,fNEtaBins};
 343         Double_t minHisto[4] = {-TMath::Pi()/2.,-2*fMaxPlotEta};
 344         Double_t maxHisto[4] = {3.*TMath::Pi()/2,2*fMaxPlotEta};
 345
 346     // --- EVENT QA PLOTS ---
 347
 348         fCentrality = new TH1F("fCentrality","Centrality;Centrality;N_{evt}",100,0,100);
 349         fHistoList->Add(fCentrality);
 350         fVertexZ = new TH1F("fVertexZ","Vertex Z position;z (cm);N_{evt}",30,-15,15);
 351         fHistoList->Add(fVertexZ);
 352
 353     // --- TRACK QA PLOTS ---
 354
 355     fDCA = new TH2F("fDCA","DCA positions TPC only cuts;xy (cm);z (cm)",100,-5,5,100,-5,5);
 356     fHistoList->Add(fDCA);
 357     fDCAZoomed = new TH2F("fDCAZoomed","DCA positions TPC only cuts;xy (cm);z (cm)",100,-.5,.5,100,-.5,.5);
 358     fHistoList->Add(fDCAZoomed);
 359     fDCAZoomedTwice = new TH2F("fDCAZoomedTwice","DCA positions TPC only cuts;xy (cm);z (cm)",100,-.05,.05,100,-.05,.05);
 360     fHistoList->Add(fDCAZoomedTwice);
 361     fDCACut = new TH2F("fDCACut","DCA positions after DCA cut;xy (cm);z (cm)",100,-5,5,100,-5,5);
 362     fHistoList->Add(fDCACut);
 363     fDCAZoomedCut = new TH2F("fDCAZoomedCut","DCA positions after DCA cut;xy (cm);z (cm)",100,-.5,.5,100,-.5,.5);
 364     fHistoList->Add(fDCAZoomedCut);
 365     fDCAZoomedTwiceCut = new TH2F("fDCAZoomedTwiceCut","DCA positions after DCA cut;xy (cm);z (cm)",100,-.05,.05,100,-.05,.05);
 366     fHistoList->Add(fDCAZoomedTwiceCut);
 367
 368     fITSHits = new TH1F("fITSHits","ITS hits",3,0,3);
 369     (fITSHits->GetXaxis())->SetBinLabel(1,"No SPD hit");
 370     (fITSHits->GetXaxis())->SetBinLabel(2,"One SPD hit");
 371     (fITSHits->GetXaxis())->SetBinLabel(3,"Two SPD hits");
 372     fHistoList->Add(fITSHits);
 373
 374     // Determine the bins used in the track cut histograms.
 375     const Int_t ncuts = 5 + fDoITSCut + fDoDCACut + fDemandNoMismatch;
 376
 377     if (!fDoITSCut) {
 378         fTrackCutLabelNumbers[3]--;
 379         fTrackCutLabelNumbers[4]--;
 380         fTrackCutLabelNumbers[5]--;
 381         fTrackCutLabelNumbers[6]--;
 382         fTrackCutLabelNumbers[7]--;
 383     }
 384
 385     if (!fDoDCACut) {
 386         fTrackCutLabelNumbers[4]--;
 387         fTrackCutLabelNumbers[5]--;
 388         fTrackCutLabelNumbers[6]--;
 389         fTrackCutLabelNumbers[7]--;
 390     }
 391
 392     fTrackCutsCount = new TH1F("fTrackCutsCount","Track Cuts;;Count",ncuts,0,ncuts);
 393     (fTrackCutsCount->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[0],"Std. TPC Only");
 394     (fTrackCutsCount->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[1],Form("#eta < %3.1f",fMaxEta));
 395     if (fDoITSCut) (fTrackCutsCount->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[2],"1-2 SPD hits");
 396     if (fDoDCACut) (fTrackCutsCount->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[3],"DCA cut");
 397     (fTrackCutsCount->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[4],"p_{T} < 5.0 GeV/c");
 398     (fTrackCutsCount->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[5],"TPCpid (p_{T} < 5.0 GeV/c)");
 399     (fTrackCutsCount->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[6],"TOFpid (p_{T} < 5.0 GeV/c)");
 400     if (fDemandNoMismatch) (fTrackCutsCount->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[7],"No Mismatch (p_{T} < 5.0 GeV/c)");
 401     fHistoList->Add(fTrackCutsCount);
 402
 403     fTrackCutsPt = new TH2F("fTrackCutsPt","Track Cuts vs p_{T};;p_{T};Count",ncuts,0,ncuts,40,0,fMaxPt);
 404     (fTrackCutsPt->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[0],"Std. TPC Only");
 405     (fTrackCutsPt->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[1],Form("#eta < %3.1f",fMaxEta));
 406     if (fDoITSCut) (fTrackCutsPt->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[2],"1-2 SPD hits");
 407     if (fDoDCACut) (fTrackCutsPt->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[3],"DCA cut");
 408     (fTrackCutsPt->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[4],"p_{T} < 5.0 GeV/c");
 409     (fTrackCutsPt->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[5],"TPCpid (p_{T} < 5.0 GeV/c)");
 410     (fTrackCutsPt->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[6],"TOFpid (p_{T} < 5.0 GeV/c)");
 411     if (fDemandNoMismatch) (fTrackCutsPt->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[7],"No Mismatch (p_{T} < 5.0 GeV/c)");
 412     fHistoList->Add(fTrackCutsPt);
 413
 414     fTrackCutsEta = new TH2F("fTrackCutsEta","Track Cuts vs #eta;;#eta;Count",ncuts,0,ncuts,4*binsHisto[1],-fMaxPlotEta,fMaxPlotEta);
 415     (fTrackCutsEta->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[0],"Std. TPC Only");
 416     (fTrackCutsEta->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[1],Form("#eta < %3.1f",fMaxEta));
 417     if (fDoITSCut) (fTrackCutsEta->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[2],"1-2 SPD hits");
 418     if (fDoDCACut) (fTrackCutsEta->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[3],"DCA cut");
 419     (fTrackCutsEta->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[4],"p_{T} < 5.0 GeV/c");
 420     (fTrackCutsEta->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[5],"TPCpid (p_{T} < 5.0 GeV/c)");
 421     (fTrackCutsEta->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[6],"TOFpid (p_{T} < 5.0 GeV/c)");
 422     if (fDemandNoMismatch) (fTrackCutsEta->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[7],"No Mismatch (p_{T} < 5.0 GeV/c)");
 423     fHistoList->Add(fTrackCutsEta);
 424
 425     fTrackCutsPhi = new TH2F("fTrackCutsPhi","Track Cuts vs #phi;;#phi;Count",ncuts,0,ncuts,72,0,2.*TMath::Pi());
 426     (fTrackCutsPhi->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[0],"Std. TPC Only");
 427     (fTrackCutsPhi->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[1],Form("#eta < %3.1f",fMaxEta));
 428     if (fDoITSCut) (fTrackCutsPhi->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[2],"1-2 SPD hits");
 429     if (fDoDCACut) (fTrackCutsPhi->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[3],"DCA cut");
 430     (fTrackCutsPhi->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[4],"p_{T} < 5.0 GeV/c");
 431     (fTrackCutsPhi->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[5],"TPCpid (p_{T} < 5.0 GeV/c)");
 432     (fTrackCutsPhi->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[6],"TOFpid (p_{T} < 5.0 GeV/c)");
 433     if (fDemandNoMismatch) (fTrackCutsPhi->GetXaxis())->SetBinLabel(fTrackCutLabelNumbers[7],"No Mismatch (p_{T} < 5.0 GeV/c)");
 434     fHistoList->Add(fTrackCutsPhi);
 435
 436     fEtaSpectrumTrig = new TH1F("fEtaSpectrumTrig","#eta Spectrum Triggers;#eta;Count",4*binsHisto[1],-fMaxPlotEta,fMaxPlotEta);
 437     fHistoList->Add(fEtaSpectrumTrig);
 438     fEtaSpectrumAssoc = new TH2F("fEtaSpectrumAssoc","#eta Spectrum Associateds;#eta;p_{T};Count",4*binsHisto[1],-fMaxPlotEta,fMaxPlotEta,10,0.,5.);
 439     fHistoList->Add(fEtaSpectrumAssoc);
 440     fPhiSpectrumAssoc = new TH2F("fPhiSpectrumAssoc","#phi Spectrum Associateds;#phi;p_{T};Count",72,0,2.*TMath::Pi(),10,0.,5.);
 441     fHistoList->Add(fPhiSpectrumAssoc);
 442
 443     // --- QA PLOTS PID ---
 444
 445         fTPCnSigmaProton = new TH2F("fTPCnSigmaProton","n#sigma plot for the TPC (Protons);p (GeV/c);n#sigma",100,0.,5.,100,-10.,10.);
 446         fHistoList->Add(fTPCnSigmaProton);
 447         fTOFnSigmaProton = new TH2F("fTOFnSigmaProton","n#sigma plot for the TOF (Protons);p (GeV/c);n#sigma",100,0.,5.,100,-10.,10.);
 448         fHistoList->Add(fTOFnSigmaProton);
 449         fTPCnSigmaPion = new TH2F("fTPCnSigmaPion","n#sigma plot for the TPC (Pions);p (GeV/c);n#sigma",100,0.,5.,100,-10.,10.);
 450         fHistoList->Add(fTPCnSigmaPion);
 451         fTOFnSigmaPion = new TH2F("fTOFnSigmaPion","n#sigma plot for the TOF (Pions);p (GeV/c);n#sigma",100,0.,5.,100,-10.,10.);
 452         fHistoList->Add(fTOFnSigmaPion);
 453         fTPCnSigmaKaon = new TH2F("fTPCnSigmaKaon","n#sigma plot for the TPC (Kaons);p (GeV/c);n#sigma",100,0.,5.,100,-10.,10.);
 454         fHistoList->Add(fTPCnSigmaKaon);
 455         fTOFnSigmaKaon = new TH2F("fTOFnSigmaKaon","n#sigma plot for the TOF (Kaons);p (GeV/c);n#sigma",100,0.,5.,100,-10.,10.);
 456         fHistoList->Add(fTOFnSigmaKaon);
 457         fTPCSignal = new TH3F("fTPCSignal","TPC Signal;#eta;p_{T} GeV/c;dE/dx",25,-fMaxEta,fMaxEta,100,0.,5.,150,0.,300.);
 458         fHistoList->Add(fTPCSignal);
 459         fTOFSignal = new TH3F("fTOFSignal","TOF Signal;#eta;p_{T} GeV/c;t",25,-fMaxEta,fMaxEta,100,0.,5.,150,10000.,25000.);
 460         fHistoList->Add(fTOFSignal);
 461
 462     // --- DI-HADRON CORRELATIONS AND MIXED EVENTS WITH TPC AND TOF SIGNALS ---
 463
 464     // Unzoomed pictures.
 465         Int_t binsTPCnoZoom[3][10] = {{100,100,100,100,100,100,100,100,100,100},
 466                                   {100,100,100,100,100,100,100,100,100,100},
 467                                   {100,100,100,100,100,100,100,100,100,100}};
 468
 469         Double_t minTPCnoZoom[3][10] = {{   -50., -50.,-30.,-30.,-30.,-30.,-30.,-30.,-30.,-30.},
 470                                     {   -100., -50.,-20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.},
 471                                     {  -200.,-150.,-50.,-30.,-20.,-20.,-20.,-20.,-20.,-20.}};
 472         Double_t maxTPCnoZoom[3][10] = {{ 150., 100., 50., 30., 25., 25., 25., 25., 25., 25.},
 473                                     {  50., 100., 40., 30., 30., 30., 30., 30., 30., 30.},
 474                                     { 100.,  50., 50., 30., 30., 30., 30., 30., 30., 30.}};
 475
 476         Int_t binsTOFnoZoom[3][10] = {{100,100,100,100,100,100,100,100,100,100},
 477                                   {100,100,100,100,100,100,100,100,100,100},
 478                                   {100,100,100,100,100,100,100,100,100,100}};
 479
 480         Double_t minTOFnoZoom[3][10] = {{-2000.,-2000.,-1000., -500., -500., -500., -500.,      -500., -500., -500.},
 481                                     {-2000.,-2000.,-2000.,-1000.,-1000.,-1000.,-1000.,-1000.,-1000.,-1000.},
 482                                     {-2000.,-2000.,-6000.,-3000.,-2000.,-1500.,-1500.,-1500.,-1500.,-1500.}};
 483         Double_t maxTOFnoZoom[3][10] = {{ 2000., 2000., 5000., 3000., 2000., 1500., 1500., 1500., 1500., 1500.},
 484                                     { 2000., 2000., 5000., 2000., 1500., 1500., 1500., 1500., 1500., 1500.},
 485                                     { 2000., 2000., 2000., 1000., 1000., 1000., 1000., 1000., 1000., 1000.}};
 486
 487     // Zoomed pictures.
 488     Int_t binsTPCZoom[3][10] = {{100,100,100,100,100,100,100,100,100,100},
 489                                 {100,100,100,100,100,100,100,100,100,100},
 490                                 {100,100,100,100,100,100,100,100,100,100}};
 491
 492         Double_t minTPCZoom[3][10] = {{   -20., -20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.},
 493                                   {      -20., -20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.},
 494                                   {   -40., -20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.,-20.}};
 495         Double_t maxTPCZoom[3][10] = {{  20.,  20., 20., 20., 20., 20., 20., 20., 20., 20.},
 496                                   {  20.,  20., 20., 20., 20., 20., 20., 20., 20., 20.},
 497                                   {  40.,  20., 20., 20., 20., 30., 30., 30., 30., 30.}};
 498
 499         Int_t binsTOFZoom[3][10] = {{100,100,100,100,100,100,100,100,100,100},
 500                                 {100,100,100,100,100,100,100,100,100,100},
 501                                 {100,100,100,100,100,100,100,100,100,100}};
 502
 503         Double_t minTOFZoom[3][10] = {{-1000.,-1000., -500., -500., -500., -400., -400., -400., -400., -400.},
 504                                   { -800., -800., -800., -800., -800., -600., -500., -500., -400., -400.},
 505                                   {-1000.,-1000.,-1000.,-1000., -800.,-1000.,-1000., -800., -700., -700.}};
 506         Double_t maxTOFZoom[3][10] = {{ 1000., 1000., 1000., 1000., 1000., 1000., 1000.,  900.,  800.,  700.},
 507                                   { 1000., 1000.,  500.,  500.,  500.,  900.,  700.,  700.,  600.,  500.},
 508                                   { 1000., 1000., 1000.,  500.,  500.,  500.,  500.,  500.,  500.,  500.}};
 509
 510         const char* species[3] = {"Pion","Kaon","Proton"};
 511
 512         for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
 513
 514                 for (Int_t iPtBin = 0; iPtBin < 10; iPtBin++) {
 515
 516             if (fZoomed) {
 517                 binsHisto[2] = binsTPCZoom[iSpecies][iPtBin];
 518                 binsHisto[3] = binsTOFZoom[iSpecies][iPtBin];
 519                 minHisto[2] = minTPCZoom[iSpecies][iPtBin];
 520                 minHisto[3] = minTOFZoom[iSpecies][iPtBin];
 521                 maxHisto[2] = maxTPCZoom[iSpecies][iPtBin];
 522                 maxHisto[3] = maxTOFZoom[iSpecies][iPtBin];
 523             } else {
 524                 binsHisto[2] = binsTPCnoZoom[iSpecies][iPtBin];
 525                 binsHisto[3] = binsTOFnoZoom[iSpecies][iPtBin];
 526                 minHisto[2] = minTPCnoZoom[iSpecies][iPtBin];
 527                 minHisto[3] = minTOFnoZoom[iSpecies][iPtBin];
 528                 maxHisto[2] = maxTPCnoZoom[iSpecies][iPtBin];
 529                 maxHisto[3] = maxTOFnoZoom[iSpecies][iPtBin];
 530             }
 531
 532             // Di-Hadron Correlations
 533             fDiHadronTPCTOF[iSpecies][iPtBin] = new THnSparseF(Form("fDiHadronTPCTOF_%s_%i",species[iSpecies],iPtBin),
 534                                                                                                                            Form("Di-Hadron Correlation (%s) %3.1f < p_{T} < %3.1f;#Delta#phi;#Delta#eta;dE/dx;t (ms);",species[iSpecies],iPtBin*0.5,(iPtBin+1)*0.5),
 535                                                                                                                            4,binsHisto,minHisto,maxHisto);
 536             fHistoList->Add(fDiHadronTPCTOF[iSpecies][iPtBin]);
 537
 538             // Mixed Events.
 539                         fMixedEventsTPCTOF[iSpecies][iPtBin] = new THnSparseF(Form("fMixedEventsTPCTOF_%s_%i",species[iSpecies],iPtBin),
 540                                                                                                                                   Form("Mixed Events (%s) %3.1f < p_{T} < %3.1f;#Delta#phi;#Delta#eta;dE/dx;t (ms);",species[iSpecies],iPtBin*0.5,(iPtBin+1)*0.5),
 541                                                                                                                                   4,binsHisto,minHisto,maxHisto);
 542                         //fHistoList->Add(fMixedEventsTPCTOF[iSpecies][iPtBin]);
 543
 544             // Inclusive spectra.
 545             fInclusiveTPCTOF[iSpecies][iPtBin] = new TH3F(Form("fInclusiveTPCTOF_%s_%i",species[iSpecies],iPtBin),
 546                                                           Form("Inclusive Spectrum (%s) %3.1f < p_{T} < %3.1f;dE/dx;t (ms);#eta",species[iSpecies],iPtBin*0.5,(iPtBin+1)*0.5),
 547                                                           binsHisto[2],minHisto[2],maxHisto[2],binsHisto[3],minHisto[3],maxHisto[3],fNEtaBins,-fMaxEta,fMaxEta);
 548             fHistoList->Add(fInclusiveTPCTOF[iSpecies][iPtBin]);
 549
 550             // Inclusive spectra with additional rapidity cut.
 551             fInclusiveTPCTOFRap[iSpecies][iPtBin] = new TH3F(Form("fInclusiveTPCTOFRap_%s_%i",species[iSpecies],iPtBin),
 552                                                              Form("Inclusive Spectrum (%s) %3.1f < p_{T} < %3.1f |Y| < %3.1f;dE/dx;t (ms);Y",species[iSpecies],iPtBin*0.5,(iPtBin+1)*0.5,fMaxRap),
 553                                                              binsHisto[2],minHisto[2],maxHisto[2],binsHisto[3],minHisto[3],maxHisto[3],20,-fMaxRap,fMaxRap);
 554             fHistoList->Add(fInclusiveTPCTOFRap[iSpecies][iPtBin]);
 555
 556
 557                 }
 558         }
 559
 560         // Correlations without PID signals
 561         fDiHadron = new TH3F("fDiHadron","Di-Hadron Correlations;#Delta#phi;#Delta#eta;p_{T_assoc}",binsHisto[0],minHisto[0],maxHisto[0],binsHisto[1],minHisto[1],maxHisto[1],10,0.,5.);
 562     fDiHadron->Sumw2();
 563         fHistoList->Add(fDiHadron);
 564
 565         fMixedEvents = new TH3F("fMixedEvents","Mixed Events;#Delta#phi;#Delta#eta;p_{T_assoc}",binsHisto[0],minHisto[0],maxHisto[0],binsHisto[1],minHisto[1],maxHisto[1],10,0.,5.);
 566         fMixedEvents->Sumw2();
 567         fHistoList->Add(fMixedEvents);
 568
 569         const char* MCSpeciesName[6] = {"piplus","pimin","Kplus","Kmin","p","pbar"};
 570         const char* MCSpeciesTitle[6] = {"#pi^{+}","#pi^{-}","K^{+}","K^{-}","p","#bar{p}"};
 571         //const char* Cuts[3] = {"nocut","DCAcut","PIDandDCAcut"};
 572
 573         // Efficiency Plots (Monte Carlo)
 574         if (fMC) {
 575
 576                 for (Int_t iSpecies=0; iSpecies<6; iSpecies++) {
 577
 578             // Without rapidity cut.
 579             fPtEtaDistrDataPrim[iSpecies]=new TH2F(Form("fPtEtaDistrDataPrim_%s",MCSpeciesName[iSpecies]),
 580                                                        Form("p_{T}-#eta distribution Data Primaries (%s);p_{T};#eta;Count",MCSpeciesTitle[iSpecies]),
 581                                                        20,0.0,5.0,fNEtaBins,-fMaxEta,fMaxEta);
 582             fPtEtaDistrDataPrim[iSpecies]->Sumw2();
 583                         fHistoList->Add(fPtEtaDistrDataPrim[iSpecies]);
 584
 585             fPtEtaDistrDataSec[iSpecies]=new TH2F(Form("fPtEtaDistrDataSec_%s",MCSpeciesName[iSpecies]),
 586                                                       Form("p_{T}-#eta distribution Data Secondaries (%s);p_{T};#eta;Count",MCSpeciesTitle[iSpecies]),
 587                                                       20,0.0,5.0,fNEtaBins,-fMaxEta,fMaxEta);
 588             fPtEtaDistrDataSec[iSpecies]->Sumw2();
 589                         fHistoList->Add(fPtEtaDistrDataSec[iSpecies]);
 590
 591                         fPtEtaDistrMCPrim[iSpecies]=new TH2F(Form("fPtEtaDistrMCPrim_%s",MCSpeciesName[iSpecies]),
 592                                                                                           Form("p_{T}-#eta distribution MC Primaries (%s);p_{T};#eta;Count",MCSpeciesTitle[iSpecies]),
 593                                                                                           20,0.0,5.0,fNEtaBins,-fMaxEta,fMaxEta);
 594                         fPtEtaDistrMCPrim[iSpecies]->Sumw2();
 595                         fHistoList->Add(fPtEtaDistrMCPrim[iSpecies]);
 596
 597                         fPtEtaDistrMCSec[iSpecies]=new TH2F(Form("fPtEtaDistrMCSec_%s",MCSpeciesName[iSpecies]),
 598                                                                                          Form("p_{T}-#eta distribution MC Secondaries (%s);p_{T};#eta;Count",MCSpeciesTitle[iSpecies]),
 599                                                                                          20,0.0,5.0,fNEtaBins,-fMaxEta,fMaxEta);
 600                         fPtEtaDistrMCSec[iSpecies]->Sumw2();
 601                         fHistoList->Add(fPtEtaDistrMCSec[iSpecies]);
 602
 603                         //fDiHadronMC[iSpecies]=new TH3F(Form("fDiHadronMC_%s",MCSpeciesName[iSpecies]),
 604                         //                                                         Form("Di-Hadron Correlations MC (%s);#Delta#phi;#Delta#eta;p_{T_assoc}",MCSpeciesTitle[iSpecies]),
 605                         //                                                         binsHisto[0],minHisto[0],maxHisto[0],binsHisto[1],minHisto[1],maxHisto[1],10,0.,5.);
 606                         //fHistoList->Add(fDiHadronMC[iSpecies]);
 607
 608             // With rapidity cut.
 609             fPtRapDistrDataPrimRapCut[iSpecies]=new TH2F(Form("fPtRapDistrDataPrimRapCut_%s",MCSpeciesName[iSpecies]),
 610                                                           Form("p_{T}-Y distribution Data Primaries (%s) |Y| < %3.1f;p_{T};Y;Count",MCSpeciesTitle[iSpecies],fMaxRap),
 611                                                           20,0.0,5.0,20,-fMaxRap,fMaxRap);
 612             fPtRapDistrDataPrimRapCut[iSpecies]->Sumw2();
 613                         fHistoList->Add(fPtRapDistrDataPrimRapCut[iSpecies]);
 614
 615             fPtRapDistrDataSecRapCut[iSpecies]=new TH2F(Form("fPtRapDistrDataSecRapCut_%s",MCSpeciesName[iSpecies]),
 616                                                       Form("p_{T}-Y distribution Data Secondaries (%s) |Y| < %3.1f;p_{T};Y;Count",MCSpeciesTitle[iSpecies],fMaxRap),
 617                                                       20,0.0,5.0,20,-fMaxRap,fMaxRap);
 618             fPtRapDistrDataSecRapCut[iSpecies]->Sumw2();
 619                         fHistoList->Add(fPtRapDistrDataSecRapCut[iSpecies]);
 620
 621                         fPtRapDistrMCPrimRapCut[iSpecies]=new TH2F(Form("fPtRapDistrMCPrimRapCut_%s",MCSpeciesName[iSpecies]),
 622                                                                                           Form("p_{T}-Y distribution MC Primaries (%s) |Y| < %3.1f;p_{T};Y;Count",MCSpeciesTitle[iSpecies],fMaxRap),
 623                                                                                           20,0.0,5.0,20,-fMaxRap,fMaxRap);
 624                         fPtRapDistrMCPrimRapCut[iSpecies]->Sumw2();
 625                         fHistoList->Add(fPtRapDistrMCPrimRapCut[iSpecies]);
 626
 627                         fPtRapDistrMCSecRapCut[iSpecies]=new TH2F(Form("fPtRapDistrMCSecRapCut_%s",MCSpeciesName[iSpecies]),
 628                                                                                          Form("p_{T}-Y distribution MC Secondaries (%s) |Y| < %3.1f;p_{T};Y;Count",MCSpeciesTitle[iSpecies],fMaxRap),
 629                                                                                          20,0.0,5.0,20,-fMaxRap,fMaxRap);
 630                         fPtRapDistrMCSecRapCut[iSpecies]->Sumw2();
 631                         fHistoList->Add(fPtRapDistrMCSecRapCut[iSpecies]);
 632
 633                 }
 634         }
 635
 636         PostData(1, fHistoList);
 637
 638 }
 639
 640 //_____________________________________________________________________________
 641 Int_t AliAnalysisTaskDiHadronPID::ClassifyTrack(AliAODTrack *track)
 642
 643 {
 644         //
 645     // This function both classifies tracks, and fills the track QA histograms.
 646     //
 647     // Classifies the track in:
 648     //  0 -> Not Useful,
 649     //  1 -> Associated,
 650     //  2 -> Trigger.
 651     //
 652     // IDEA: later we can do this with filterbits.
 653     //
 654     // The following track cuts are applied for triggers:
 655     //
 656     // 1a) pT > 5.0 GeV/c, pT < fPtMax,
 657     //  2) StandardTPCOnlyTrackCuts (filterbit 7 for AOD's),
 658     //  3) eta < fMaxEta,
 659     //  4) ITS track cut, a track is only selected if it has at least one hit in the SPD,
 660     //     that is, in one of the first two layers of the ITS. (can be switched on and off
 661     //     using (Bool_t)fDoITSCut),
 662     //  5) DCA track cut, of all tracks with at least one SPD hit, the DCA is constrained as
 663     //     a function of pT, in order to remove lots of the secondaries. (can be switched on and
 664     //     off by using (Bool_t)fDoDCACut),
 665     //
 666     // For associateds all the same track cuts are applied, except 1a is
 667     // replaced by 2b. Moreover the following cuts are applied too:
 668     //
 669     // 1b) pT < 5.0 GeV/c,
 670     //  6) TPCpid,
 671     //  7) TOFpid,
 672     //  8) no TPCTOF mismatch.
 673     //
 674
 675     // Check if a track is supplied.
 676     if (!track) {
 677         if (fVerbose>0) cout<<"AliAnalysisTaskDiHadronPID::SelectTrack -> ERROR: No track found."<<endl;
 678         return 0;
 679     }
 680
 681     // Get basic track information.
 682     Int_t classification=0;
 683     Double_t pt = track->Pt();
 684     Double_t eta = track->Eta();
 685     Double_t phi = track->Phi();
 686
 687     // 1) pT cut: First separation between triggers and associateds.
 688     if (pt<5.0) classification=1;
 689     if ((pt>5.0)&&(pt<fMaxPt)) classification = 2;
 690     if (!classification) return 0;
 691
 692     // 2) StandardTPCOnlyTrackCuts.
 693     if (!(track->TestFilterMask(1<<7))) {
 694                 //if (fVerbose>3) cout<<"Track Ignored: Did Not pass filterbit."<<endl;
 695                 return 0;
 696         }
 697
 698         if (fVerbose>3) {
 699                 cout<<endl;
 700                 cout<<"pt: "<<pt<<" eta: "<<eta<<" phi: "<<phi<<endl;
 701         }
 702
 703     fTrackCutsCount->Fill(fTrackCutLabelNumbers[0]-.5);
 704     fTrackCutsPt->Fill(fTrackCutLabelNumbers[0]-.5,pt);
 705     fTrackCutsEta->Fill(fTrackCutLabelNumbers[0]-.5,eta);
 706     fTrackCutsPhi->Fill(fTrackCutLabelNumbers[0]-.5,phi);
 707
 708     // 3) eta cut.
 709     if (TMath::Abs(eta)>fMaxEta) {
 710                 if (fVerbose>3) cout<<"Track Ignored: Eta too large."<<endl;
 711                 return 0;
 712     }
 713
 714     fTrackCutsCount->Fill(fTrackCutLabelNumbers[1]-.5);
 715     fTrackCutsPt->Fill(fTrackCutLabelNumbers[1]-.5,pt);
 716     fTrackCutsEta->Fill(fTrackCutLabelNumbers[1]-.5,eta);
 717     fTrackCutsPhi->Fill(fTrackCutLabelNumbers[1]-.5,phi);
 718
 719     // Obtaining ITS information.
 720     AliAODTrack* globaltrack = GetGlobalTrack(track);
 721     if (!globaltrack) {
 722       return 0;
 723     }
 724
 725     Bool_t ITSLayerHit[6];
 726     for (Int_t iITSLayer=0; iITSLayer<6; iITSLayer++) {
 727         ITSLayerHit[iITSLayer] = globaltrack->HasPointOnITSLayer(iITSLayer);
 728     }
 729     Int_t SPDHits=ITSLayerHit[0]+ITSLayerHit[1];
 730
 731         if (fVerbose>3) cout<<"SPD hits: "<<SPDHits<<endl;
 732
 733     // Fill the ITS hist.
 734     fITSHits->Fill(SPDHits+0.5);
 735
 736     // 4) ITS cut.
 737     if (fDoITSCut) {
 738         if (!SPDHits) {
 739                         if (fVerbose>3) cout<<"Track Ignored: Not enough SPD hits."<<endl;
 740                         return 0;
 741                 }
 742         fTrackCutsCount->Fill(fTrackCutLabelNumbers[2]-.5);
 743         fTrackCutsPt->Fill(fTrackCutLabelNumbers[2]-.5,pt);
 744         fTrackCutsEta->Fill(fTrackCutLabelNumbers[2]-.5,eta);
 745         fTrackCutsPhi->Fill(fTrackCutLabelNumbers[2]-.5,phi);
 746     }
 747
 748     // Propagate the global track to the DCA.
 749     Double_t PosAtDCA[2] = {-999,-999};
 750     Double_t covar[3] = {-999,-999,-999};
 751     globaltrack->PropagateToDCA(fAODVertex,fAODEvent->GetMagneticField(),100.,PosAtDCA,covar);
 752
 753     // Fill the DCA hist (before cut)
 754     fDCA->Fill(PosAtDCA[0],PosAtDCA[1]);
 755     fDCAZoomed->Fill(PosAtDCA[0],PosAtDCA[1]);
 756     fDCAZoomedTwice->Fill(PosAtDCA[0],PosAtDCA[1]);
 757
 758     // 5) DCA cut (See R_AA paper).
 759     Double_t DCAcutvalue[2];
 760     DCAcutvalue[0] = 0.018 + 0.035*TMath::Power(pt,-1.01);
 761     DCAcutvalue[1] = 2.;
 762     if (fDoDCACut) {
 763         if (SPDHits&&((TMath::Abs(PosAtDCA[0])>DCAcutvalue[0])||(TMath::Abs(PosAtDCA[1])>DCAcutvalue[1]))) {
 764                         if (fVerbose>3) cout<<"Track Ignored: Enough SPD hits, but out of DCA range."<<endl;
 765                         return 0;
 766                 }
 767         fTrackCutsCount->Fill(fTrackCutLabelNumbers[3]-.5);
 768         fTrackCutsPt->Fill(fTrackCutLabelNumbers[3]-.5,pt);
 769         fTrackCutsEta->Fill(fTrackCutLabelNumbers[3]-.5,eta);
 770         fTrackCutsPhi->Fill(fTrackCutLabelNumbers[3]-.5,phi);
 771
 772         // Fill the DCA hist (after cut)
 773         fDCACut->Fill(PosAtDCA[0],PosAtDCA[1]);
 774         fDCAZoomedCut->Fill(PosAtDCA[0],PosAtDCA[1]);
 775         fDCAZoomedTwiceCut->Fill(PosAtDCA[0],PosAtDCA[1]);
 776     }
 777
 778     // Now all the common cuts have been performed. Tracks identified as a trigger will
 779     // be returned. Note that they will not appear in the histograms after track cut 5.
 780     if (classification==2) return 2;
 781
 782     // Now we're left with only tracks with pT < 5.0 GeV/c.
 783     fTrackCutsCount->Fill(fTrackCutLabelNumbers[4]-.5);
 784     fTrackCutsPt->Fill(fTrackCutLabelNumbers[4]-.5,pt);
 785     fTrackCutsEta->Fill(fTrackCutLabelNumbers[4]-.5,eta);
 786     fTrackCutsPhi->Fill(fTrackCutLabelNumbers[4]-.5,phi);
 787
 788     // Obtain the status of the track.
 789     ULong_t status = globaltrack->GetStatus();
 790
 791     // 6) TPCpid
 792     if (!((status&AliAODTrack::kTPCpid)==AliAODTrack::kTPCpid)) {
 793                 if (fVerbose>3) cout<<"Track Ignored: No TPC pid."<<endl;
 794                 return 0;
 795         }
 796
 797     fTrackCutsCount->Fill(fTrackCutLabelNumbers[5]-.5);
 798     fTrackCutsPt->Fill(fTrackCutLabelNumbers[5]-.5,pt);
 799     fTrackCutsEta->Fill(fTrackCutLabelNumbers[5]-.5,eta);
 800     fTrackCutsPhi->Fill(fTrackCutLabelNumbers[5]-.5,phi);
 801
 802     // 7) TOFpid
 803     if (!((status&AliAODTrack::kTOFpid)==AliAODTrack::kTOFpid)) {
 804                 if (fVerbose>3) cout<<"Track Ignored: No TOF pid."<<endl;
 805                 return 0;
 806         }
 807
 808     fTrackCutsCount->Fill(fTrackCutLabelNumbers[6]-.5);
 809     fTrackCutsPt->Fill(fTrackCutLabelNumbers[6]-.5,pt);
 810     fTrackCutsEta->Fill(fTrackCutLabelNumbers[6]-.5,eta);
 811     fTrackCutsPhi->Fill(fTrackCutLabelNumbers[6]-.5,phi);
 812
 813     // 8) TPC, TOF mismatch.
 814     if (fDemandNoMismatch) {
 815         if ((status&AliAODTrack::kTOFmismatch)==AliAODTrack::kTOFmismatch) {
 816                         if (fVerbose>3) cout<<"Track Ignored: TOF mismatch found."<<endl;
 817                         return 0;
 818                 }
 819         fTrackCutsCount->Fill(fTrackCutLabelNumbers[7]-.5);
 820         fTrackCutsPt->Fill(fTrackCutLabelNumbers[7]-.5,pt);
 821         fTrackCutsEta->Fill(fTrackCutLabelNumbers[7]-.5,eta);
 822         fTrackCutsPhi->Fill(fTrackCutLabelNumbers[7]-.5,phi);
 823     }
 824
 825     // All tracks which made it up to here are classified as associateds.
 826
 827     // Return associated.
 828     return 1;
 829
 830 }
 831
 832 //____________________________________________________________________________
 833 Bool_t AliAnalysisTaskDiHadronPID::SelectEvent(AliAODVertex* vertex)
 834
 835 {
 836
 837         //
 838         // Event Selection.
 839         //
 840
 841         Double_t primVtx[3];
 842         vertex->GetXYZ(primVtx);
 843         if (TMath::Sqrt(primVtx[0]*primVtx[0] + primVtx[1]*primVtx[1])>1. || TMath::Abs(primVtx[2])>10.) {
 844                 if (fVerbose>2) cout << "AliAnalysisTaskDiHadronPID::SelectEvent -> Vertex Out of Range." << endl;
 845         if (fVerbose>2) cout << "AliAnalysisTaskDiHadronPID::SelectEvent -> Event not selected." << endl;
 846                 return kFALSE;
 847         }
 848         if (fVerbose>2) cout << "AliAnalysisTaskDiHadronPID::SelectEvent -> Vertex is OK." << endl;
 849
 850     // We also wish to make centrality cut, but only if run on Pb+Pb.
 851     if (fBeamType=="PbPb") {
 852         Double_t cent = fAODHeader->GetCentrality();
 853         if (cent>fCentralityCutMin||cent<fCentralityCutMax) {
 854             if (fVerbose>2) cout<<"AliAnalysisTaskDiHadronPID::SelectEvent -> Event did not pass centaltity cut."<<endl;
 855             return kFALSE;
 856         }
 857     }
 858
 859     if (fVerbose>2) cout << "AliAnalysisTaskDiHadronPID::SelectEvent -> Event selected." << endl;
 860         return kTRUE;
 861
 862 }
 863
 864 //_____________________________________________________________________________
 865 void AliAnalysisTaskDiHadronPID::FillGlobalTracksArray() {
 866
 867         // Initialize the mapping for corresponding PID tracks. (see
 868         // GetGlobalTrack(AliAODTrack* track)).
 869         //
 870
 871         if (!fAODEvent) {
 872         if (fVerbose>0) cout << "AliAnalysisTaskDiHadronPID::FillGlobalTracksArray -> ERROR: fAODEvent not set." << endl;
 873                 return;
 874         }
 875
 876         fGlobalTracks = new TObjArray();
 877         AliAODTrack* track = 0x0;
 878
 879         for (Int_t iTrack = 0; iTrack < fAODEvent->GetNumberOfTracks(); iTrack++) {
 880
 881                 track = fAODEvent->GetTrack(iTrack);
 882
 883                 // I.e., if it does NOT pass the filtermask.
 884                 if (!(track->TestFilterMask(1<<7))) {
 885             if (track->GetID()>-1) fGlobalTracks->AddAtAndExpand(track,track->GetID());
 886             //if (track->GetID()<1) cout<<"Track ID: "<<track->GetID()<<" Partner ID: "<<(-track->GetID()-1)<<endl;
 887                 }
 888
 889         }
 890
 891 }
 892
 893 //_____________________________________________________________________________
 894 AliAODTrack* AliAnalysisTaskDiHadronPID::GetGlobalTrack(AliAODTrack* track) {
 895
 896         //
 897         // Returns the "parner track" of track, which contains the pid information.
 898         //
 899
 900         AliAODTrack* partner = 0x0;
 901
 902     partner = (AliAODTrack*)(fGlobalTracks->At(-track->GetID()-1));
 903
 904         if (!partner&&(fVerbose>3)) cout<<"AliAnalysisTaskDiHadronPID::GetGlobalTrack -> No Global Track Found!"<<endl;
 905
 906         return partner;
 907
 908 }
 909
 910 //_____________________________________________________________________________
 911 Double_t AliAnalysisTaskDiHadronPID::PhiRange(Double_t DPhi)
 912
 913 {
 914         //
 915         // Puts the argument in the range [-pi/2,3 pi/2].
 916         //
 917
 918         if (DPhi < -TMath::Pi()/2) DPhi += 2*TMath::Pi();
 919         if (DPhi > 3*TMath::Pi()/2) DPhi -= 2*TMath::Pi();
 920
 921         return DPhi;
 922
 923 }
 924
 925 //_____________________________________________________________________________
 926 Int_t AliAnalysisTaskDiHadronPID::ConvertPdgCode(Int_t pdgcode) {
 927
 928     //
 929     // Converts the pdg code to the bin number (0-5)
 930     //
 931
 932     if (pdgcode==211) return 0; // Pi +
 933     if (pdgcode==-211) return 1; // Pi -
 934     if (pdgcode==321) return 2; // K +
 935     if (pdgcode==-321) return 3; // K -
 936     if (pdgcode==2212) return 4; // p +
 937     if (pdgcode==-2212) return 5; // p -
 938
 939     return -999; // Any other particle.
 940
 941 }
 942
 943 //_____________________________________________________________________________
 944 void AliAnalysisTaskDiHadronPID::UserExec(Option_t *)
 945
 946 {
 947         //
 948         // UserExec.
 949         //
 950
 951         // Input the event.
 952         fAODEvent = dynamic_cast<AliAODEvent*>(InputEvent());
 953         if (!fAODEvent) {
 954                 if (fVerbose>0) cout << "AliAnalysisTaskDiHadronPID::UserExec -> ERROR: No AliAODEvent pointer could be created." << endl;
 955                 return;
 956         }
 957
 958         // Get the event header.
 959         fAODHeader = fAODEvent->GetHeader();
 960         if (!fAODHeader) {
 961                 if (fVerbose>0) cout << "AliAnalysisTaskDiHadronPID::UserExec -> ERROR: No AliAODHeader pointer could be created."<<endl;
 962                 return;
 963         }
 964
 965         // Get the event vertex.
 966         fAODVertex = fAODEvent->GetPrimaryVertex();
 967         if (!fAODVertex) {
 968                 if (fVerbose>0) cout << "AliAnalysisTaskDiHadronPID::UserExec -> ERROR: No AliAODVertex pointer could be created." << endl;
 969                 return;
 970         }
 971
 972     // Get the MC tracks if ran on MC.
 973     if (fMC) {
 974         fMCTracks = dynamic_cast<TClonesArray*>(fAODEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 975         if (!fMCTracks) {
 976             if (fVerbose>0) {
 977                 cout<<"AliAnalysisTaskDiHadronPID::UserExec -> ERROR: No MC particles found."<<endl;
 978                 return;
 979             }
 980         }
 981     }
 982
 983         AliTPCPIDResponse& TPCPIDResponse = fPIDResponse->GetTPCResponse();
 984
 985     // See if the event passes the event selection.
 986     if (!SelectEvent(fAODVertex)) return;
 987
 988         // Display basic event information.
 989         if ((fVerbose>2)) cout << endl;
 990         if ((fVerbose>2)&&(fBeamType=="PbPb")) cout << "Event centrality: " << fAODHeader->GetCentrality() << endl;
 991         if ((fVerbose>2)) cout << "Event Vertex Z: " << fAODVertex->GetZ() << endl;
 992         if ((fVerbose>2)) cout << "Event tracks in AOD: " << fAODEvent->GetNumberOfTracks() << endl;
 993         if ((fVerbose>2)) cout << endl;
 994
 995     // Filling Event QA plots.
 996         if (fBeamType=="PbPb") fCentrality->Fill(fAODHeader->GetCentrality());
 997         fVertexZ->Fill(fAODVertex->GetZ());
 998
 999         // Fill the TObjArray which holds Global tracks.
1000         FillGlobalTracksArray();
1001
1002         // Create object arrays for triggers and associateds.
1003         TObjArray *triggers             = new TObjArray();
1004         TObjArray *associateds  = new TObjArray();
1005
1006     // Loop over MC tracks to fill the spectra if ran on MC.
1007     if (fMC) {
1008         for (Int_t iTrack=0; iTrack<fMCTracks->GetEntriesFast(); iTrack++) {
1009
1010             AliAODMCParticle* mcparticle = (AliAODMCParticle*)fMCTracks->At(iTrack);
1011             Double_t mcPt = mcparticle->Pt();
1012             Double_t mcEta = mcparticle->Eta();
1013             //Double_t mcPhi = mcparticle->Phi();
1014             Double_t mcY = mcparticle->Y();
1015             Int_t mcPDG = mcparticle->PdgCode();
1016             Int_t mcSpeciesBin = ConvertPdgCode(mcPDG);
1017             if (mcSpeciesBin==-999) continue;
1018
1019             if (mcparticle->IsPhysicalPrimary()) {
1020                 fPtEtaDistrMCPrim[mcSpeciesBin]->Fill(mcPt,mcEta);
1021                 if (TMath::Abs(mcY)<fMaxRap) fPtRapDistrMCPrimRapCut[mcSpeciesBin]->Fill(mcPt,mcY);
1022             } else {
1023                 fPtEtaDistrMCSec[mcSpeciesBin]->Fill(mcPt,mcEta);
1024                 if (TMath::Abs(mcY)<fMaxRap) fPtRapDistrMCSecRapCut[mcSpeciesBin]->Fill(mcPt,mcY);
1025             }
1026         }
1027     }
1028
1029     // In this loop the triggers and associateds will be identified, track QA and PID QA histograms will be filled.
1030         for (Int_t iTrack = 0; iTrack < fAODEvent->GetNumberOfTracks(); iTrack++) {
1031
1032         // Obtain a pointer to the track.
1033                 fAODTrack = fAODEvent->GetTrack(iTrack);
1034         if (!fAODTrack&&(fVerbose>0)) {
1035             cout << "AliAnalysisTaskDiHadronPID::UserExec -> ERROR: Track object not found." << endl;
1036             continue;
1037         }
1038
1039                 //if (TMath::Abs(fAODTrack->Eta())>0.8&&fAODTrack->Pt()>0.5&&TMath::Abs(fAODTrack->Y(AliAODTrack::kProton))<0.5) cout<<"Eta: "<<fAODTrack->Eta()<<" Pt: "<<fAODTrack->Pt()<<" Ypion: "<<fAODTrack->Y(AliAODTrack::kPion)<<" Ykaon: "<<fAODTrack->Y(AliAODTrack::kKaon)<<" Yproton: "<<fAODTrack->Y(AliAODTrack::kProton)<<endl;
1040
1041         // Find the track classification.
1042         Int_t tracktype = ClassifyTrack(fAODTrack);
1043
1044         if (tracktype==0) {
1045                         continue;
1046                 }
1047
1048         if (tracktype==1) {
1049                         if (fVerbose>3) cout<<"Track added to associated buffer."<<endl;
1050             associateds->AddLast(fAODTrack);
1051             fEtaSpectrumAssoc->Fill(fAODTrack->Eta(),fAODTrack->Pt());
1052             fPhiSpectrumAssoc->Fill(fAODTrack->Phi(),fAODTrack->Pt());
1053         }
1054
1055         if (tracktype==2) {
1056                         if (fVerbose>3) cout<<"Track added to trigger buffer."<<endl;
1057             triggers->AddLast(fAODTrack);
1058             fEtaSpectrumTrig->Fill(fAODTrack->Eta());
1059
1060         }
1061
1062         if (tracktype==1) {
1063
1064             // Fill the PID QA histograms for the associateds
1065             AliAODTrack* globaltrack = GetGlobalTrack(fAODTrack);
1066             Double_t mom, nSigma;
1067
1068                         Double_t pt = fAODTrack->Pt();
1069                         Double_t eta = fAODTrack->Eta();
1070                         const Int_t ptbin = (Int_t)(2*pt);
1071
1072                         //cout<<pt<<" "<<ptbin<<endl;
1073
1074             mom = globaltrack->GetTPCmomentum();
1075             nSigma = fPIDResponse->NumberOfSigmasTPC(globaltrack,AliPID::kProton);
1076             fTPCnSigmaProton->Fill(mom,nSigma);
1077             nSigma = fPIDResponse->NumberOfSigmasTPC(globaltrack,AliPID::kPion);
1078             fTPCnSigmaPion->Fill(mom,nSigma);
1079             nSigma = fPIDResponse->NumberOfSigmasTPC(globaltrack,AliPID::kKaon);
1080             fTPCnSigmaKaon->Fill(mom,nSigma);
1081
1082             fTPCSignal->Fill(eta,pt,globaltrack->GetTPCsignal());
1083
1084             mom =globaltrack->P();
1085             nSigma = fPIDResponse->NumberOfSigmasTOF(globaltrack,AliPID::kProton);
1086             fTOFnSigmaProton->Fill(mom,nSigma);
1087             nSigma = fPIDResponse->NumberOfSigmasTOF(globaltrack,AliPID::kPion);
1088             fTOFnSigmaPion->Fill(mom,nSigma);
1089             nSigma = fPIDResponse->NumberOfSigmasTOF(globaltrack,AliPID::kKaon);
1090             fTOFnSigmaKaon->Fill(mom,nSigma);
1091
1092             fTOFSignal->Fill(eta,pt,globaltrack->GetTOFsignal());
1093
1094                         // Fill the inclusives.
1095                         Double_t TPCmom = globaltrack->GetTPCmomentum();
1096                         Double_t TPCsignal = globaltrack->GetTPCsignal();
1097                         Double_t expectedTPCsignalPion = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kPion);
1098                         Double_t expectedTPCsignalKaon = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kKaon);
1099                         Double_t expectedTPCsignalProton = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kProton);
1100
1101                         Double_t TOFsignal = globaltrack->GetTOFsignal();
1102                         Double_t times[AliPID::kSPECIES];
1103                         globaltrack->GetIntegratedTimes(times);
1104                         Double_t expectedTOFsignalPion = times[AliPID::kPion];
1105                         Double_t expectedTOFsignalKaon = times[AliPID::kKaon];
1106                         Double_t expectedTOFsignalProton = times[AliPID::kProton];
1107
1108                         Double_t TPCsignalSubtracted = TPCsignal - expectedTPCsignalPion;
1109                         Double_t TOFsignalSubtracted = TOFsignal - expectedTOFsignalPion;
1110                         fInclusiveTPCTOF[0][ptbin]->Fill(TPCsignalSubtracted,TOFsignalSubtracted,eta);
1111                         if (fAODTrack->Y(AliAODTrack::kPion)<fMaxRap) fInclusiveTPCTOFRap[0][ptbin]->Fill(TPCsignalSubtracted,TOFsignalSubtracted,fAODTrack->Y(AliAODTrack::kPion));
1112
1113                         TPCsignalSubtracted = TPCsignal - expectedTPCsignalKaon;
1114                         TOFsignalSubtracted = TOFsignal - expectedTOFsignalKaon;
1115                         fInclusiveTPCTOF[1][ptbin]->Fill(TPCsignalSubtracted,TOFsignalSubtracted,eta);
1116                         if (fAODTrack->Y(AliAODTrack::kKaon)<fMaxRap) fInclusiveTPCTOFRap[1][ptbin]->Fill(TPCsignalSubtracted,TOFsignalSubtracted,fAODTrack->Y(AliAODTrack::kKaon));
1117
1118                         TPCsignalSubtracted = TPCsignal - expectedTPCsignalProton;
1119                         TOFsignalSubtracted = TOFsignal - expectedTOFsignalProton;
1120                         fInclusiveTPCTOF[2][ptbin]->Fill(TPCsignalSubtracted,TOFsignalSubtracted,eta);
1121                         if (fAODTrack->Y(AliAODTrack::kProton)<fMaxRap) fInclusiveTPCTOFRap[2][ptbin]->Fill(TPCsignalSubtracted,TOFsignalSubtracted,fAODTrack->Y(AliAODTrack::kProton));
1122
1123             // Fill the MC reconstructed spectra histograms.
1124             if (fMC) {
1125
1126                 Int_t aodlabel = TMath::Abs(fAODTrack->GetLabel());
1127                 AliAODMCParticle* mcparticle = (AliAODMCParticle*)fMCTracks->At(aodlabel);
1128                 Int_t mcPDG = mcparticle->PdgCode();
1129                 Int_t mcSpeciesBin = ConvertPdgCode(mcPDG);
1130
1131                 // Only fill hisotos for pions, kaons and protons.
1132                 if (mcSpeciesBin!=-999) {
1133                     Double_t dataPt = fAODTrack->Pt();
1134                     Double_t dataY=-999;
1135                     if (mcSpeciesBin==0||mcSpeciesBin==1) {dataY = fAODTrack->Y(AliAODTrack::kPion);}
1136                     if (mcSpeciesBin==2||mcSpeciesBin==3) {dataY = fAODTrack->Y(AliAODTrack::kKaon);}
1137                     if (mcSpeciesBin==4||mcSpeciesBin==5) {dataY = fAODTrack->Y(AliAODTrack::kProton);}
1138
1139                                         //if (TMath::Abs(fAODTrack->Eta())>0.8) cout<<"Eta: "<<fAODTrack->Eta()<<" Ypion: "<<fAODTrack->Y(AliAODTrack::kPion)<<" Ykaon: "<<fAODTrack->Y(AliAODTrack::kKaon)<<" Yproton: "<<fAODTrack->Y(AliAODTrack::kProton)<<endl;
1140
1141                     if (mcparticle->IsPhysicalPrimary()) {
1142                         fPtEtaDistrDataPrim[mcSpeciesBin]->Fill(dataPt,eta);
1143                         if (TMath::Abs(dataY)<fMaxRap) fPtRapDistrDataPrimRapCut[mcSpeciesBin]->Fill(dataPt,dataY);
1144                     } else {
1145                         fPtEtaDistrDataSec[mcSpeciesBin]->Fill(dataPt,eta);
1146                         if (TMath::Abs(dataY)<fMaxRap) fPtRapDistrDataSecRapCut[mcSpeciesBin]->Fill(dataPt,dataY);
1147                     }
1148                 }
1149             }
1150
1151
1152
1153         }
1154     }
1155
1156     // In This Loop the di-hadron correlation will be made.
1157     Double_t histoFill[4];                  // {Dphi, Deta, TPC signal, TOF signal}
1158     AliAODTrack* currentTrigger = 0x0;
1159         AliAODTrack* currentAssociated = 0x0;
1160         AliAODTrack* currentAssociatedGlobal = 0x0;
1161
1162     for (Int_t iTrig = 0; iTrig < triggers->GetEntriesFast(); iTrig++){
1163
1164                 currentTrigger = (AliAODTrack*)(triggers->At(iTrig));
1165
1166                 for (Int_t iAssoc = 0; iAssoc < associateds->GetEntriesFast(); iAssoc++) {
1167
1168                         currentAssociated = (AliAODTrack*)(associateds->At(iAssoc));
1169                         currentAssociatedGlobal = GetGlobalTrack(currentAssociated);
1170
1171                         Double_t pt = currentAssociated->Pt();
1172                         histoFill[0] = PhiRange(currentTrigger->Phi() - currentAssociated->Phi());
1173                         histoFill[1] = currentTrigger->Eta() - currentAssociated->Eta();
1174
1175                         // Is there a caveat here when Pt = 5.00000000?
1176                         const Int_t ptbin = (Int_t)(2*pt);
1177                         // cout<<"pt: "<<pt<<" ptbin: "<<ptbin<<endl; // Works OK!
1178
1179                         if (currentAssociatedGlobal) {
1180
1181                 // Get TPC (expected) signals.
1182                                 Double_t TPCmom = currentAssociatedGlobal->GetTPCmomentum();
1183                 Double_t TPCsignal = currentAssociatedGlobal->GetTPCsignal();
1184                                 Double_t expectedTPCsignalPion = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kPion);
1185                 Double_t expectedTPCsignalKaon = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kKaon);
1186                                 Double_t expectedTPCsignalProton = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kProton);
1187
1188                 // Get TOF (expected) signals.
1189                 Double_t TOFsignal = currentAssociatedGlobal->GetTOFsignal();
1190                 Double_t times[AliPID::kSPECIES];
1191                 currentAssociatedGlobal->GetIntegratedTimes(times);
1192                 Double_t expectedTOFsignalPion = times[AliPID::kPion];
1193                 Double_t expectedTOFsignalKaon = times[AliPID::kKaon];
1194                 Double_t expectedTOFsignalProton = times[AliPID::kProton];
1195
1196                 // Fill the histograms.
1197                 histoFill[2] = TPCsignal - expectedTPCsignalPion;
1198                 histoFill[3] = TOFsignal - expectedTOFsignalPion;
1199                 fDiHadronTPCTOF[0][ptbin]->Fill(histoFill);
1200
1201                 histoFill[2] = TPCsignal - expectedTPCsignalKaon;
1202                 histoFill[3] = TOFsignal - expectedTOFsignalKaon;
1203                 fDiHadronTPCTOF[1][ptbin]->Fill(histoFill);
1204
1205                 histoFill[2] = TPCsignal - expectedTPCsignalProton;
1206                 histoFill[3] = TOFsignal - expectedTOFsignalProton;
1207                 fDiHadronTPCTOF[2][ptbin]->Fill(histoFill);
1208
1209                                 fDiHadron->Fill(histoFill[0],histoFill[1],pt);
1210
1211                         }
1212                 }
1213         }
1214
1215     // In this loop we calculate the mixed events.
1216     if (fCalculateMixedEvents) {
1217
1218         // Loop over the trigger buffer.
1219         if (fVerbose>3) cout << "AliAnalysisTaskDiHadronPID::UserExec -> Mixing the events with "<<fTrigBufferSize<<" triggers from the buffer." <<endl;
1220         if (fVerbose>3) cout << "AliAnalysisTaskDiHadronPID::UserExec -> Buffer size: "<<fTrigBufferIndex<<endl;
1221
1222         for (Int_t iTrig=0;iTrig<fTrigBufferSize;iTrig++) {
1223
1224             // Check if the trigger and the associated have a reconstructed
1225             // vertext no further than 2cm apart.
1226
1227             // fTrigBuffer[i][0] = z
1228             // fTrigBuffer[i][1] = phi
1229             // fTrigBuffer[i][2] = eta
1230             // fTrigBuffer[i][3] = p_t
1231
1232             if (TMath::Abs(fTrigBuffer[iTrig][0]-fAODVertex->GetZ())<fVertexZMixedEvents) {
1233
1234                 if (fVerbose>3) cout<<"AliAnalysisTaskDiHadronPID::UserExec -> Mixing with trigger Z: "<<fTrigBuffer[iTrig][0]<<", Pt: "<<fTrigBuffer[iTrig][3]<<endl;
1235
1236                 for (Int_t iAssoc = 0; iAssoc < associateds->GetEntriesFast(); iAssoc++) {
1237
1238                     currentAssociated = (AliAODTrack*)(associateds->At(iAssoc));
1239                     currentAssociatedGlobal = GetGlobalTrack(currentAssociated);
1240
1241                     if (currentAssociatedGlobal) {
1242
1243                                                 Double_t pt = currentAssociated->Pt();
1244                         histoFill[0] = PhiRange(fTrigBuffer[iTrig][1] - currentAssociated->Phi());
1245                         histoFill[1] = fTrigBuffer[iTrig][2] - currentAssociated->Eta();
1246
1247                                                 //const Int_t ptbin = (Int_t)(2*pt);
1248
1249                                                 // Get TPC (expected) signals.
1250                                                 Double_t TPCmom = currentAssociatedGlobal->GetTPCmomentum();
1251                                                 Double_t TPCsignal = currentAssociatedGlobal->GetTPCsignal();
1252                                                 Double_t expectedTPCsignalPion = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kPion);
1253                                                 Double_t expectedTPCsignalKaon = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kKaon);
1254                                                 Double_t expectedTPCsignalProton = TPCPIDResponse.GetExpectedSignal(TPCmom,AliPID::kProton);
1255
1256                                                 // Get TOF (expected) signals.
1257                                                 Double_t TOFsignal = currentAssociatedGlobal->GetTOFsignal();
1258                                                 Double_t times[AliPID::kSPECIES];
1259                                                 currentAssociatedGlobal->GetIntegratedTimes(times);
1260                                                 Double_t expectedTOFsignalPion = times[AliPID::kPion];
1261                                                 Double_t expectedTOFsignalKaon = times[AliPID::kKaon];
1262                                                 Double_t expectedTOFsignalProton = times[AliPID::kProton];
1263
1264                                                 // Fill the histograms.
1265                                                 histoFill[2] = TPCsignal - expectedTPCsignalPion;
1266                                                 histoFill[3] = TOFsignal - expectedTOFsignalPion;
1267                                                 //if (ptbin==1) fMixedEventsTPCTOF[0][ptbin]->Fill(histoFill);
1268
1269                                                 histoFill[2] = TPCsignal - expectedTPCsignalKaon;
1270                                                 histoFill[3] = TOFsignal - expectedTOFsignalKaon;
1271                                                 //if (ptbin==1) fMixedEventsTPCTOF[1][ptbin]->Fill(histoFill);
1272
1273                                                 histoFill[2] = TPCsignal - expectedTPCsignalProton;
1274                                                 histoFill[3] = TOFsignal - expectedTOFsignalProton;
1275                                                 //if (ptbin==1) fMixedEventsTPCTOF[2][ptbin]->Fill(histoFill);
1276
1277                                                 fMixedEvents->Fill(histoFill[0],histoFill[1],pt);
1278
1279                                                 /*
1280                                                 Double_t PionFillTPC = TPCsignal - expectedTPCsignalPion;
1281                                                 Double_t PionFillTOF = TOFsignal - expectedTOFsignalPion;
1282
1283                                                 Int_t PionMinTPC = (fDiHadronTPCTOF[0][ptbin]->GetAxis(2))->GetXmin();
1284                                                 Int_t PionMaxTPC = (fDiHadronTPCTOF[0][ptbin]->GetAxis(2))->GetXmax();
1285                                                 Int_t PionMinTOF = (fDiHadronTPCTOF[0][ptbin]->GetAxis(3))->GetXmin();
1286                                                 Int_t PionMaxTOF = (fDiHadronTPCTOF[0][ptbin]->GetAxis(3))->GetXmax();
1287
1288                                                 //cout<<PionMinTPC<<"<"<<PionFillTPC<<"<"<<PionMaxTPC<<"? "<<PionMinTOF<<"<"<<PionFillTOF<<"<"<<PionMaxTOF<<"?"<<endl;
1289                                                 if (PionFillTPC<PionMaxTPC&&PionFillTPC>PionMinTPC&&PionFillTOF<PionMaxTOF&&PionFillTOF>PionMinTOF) {
1290                                                         fMixedEventsTPCTOFCut[0]->Fill(DPhi,DEta,pt);
1291                                                         //cout<<"Yes! Histogram will be filled."<<endl;
1292                                                 }
1293
1294                                                 Double_t KaonFillTPC = TPCsignal - expectedTPCsignalKaon;
1295                                                 Double_t KaonFillTOF = TOFsignal - expectedTOFsignalKaon;
1296
1297                                                 Int_t KaonMinTPC = (fDiHadronTPCTOF[1][ptbin]->GetAxis(2))->GetXmin();
1298                                                 Int_t KaonMaxTPC = (fDiHadronTPCTOF[1][ptbin]->GetAxis(2))->GetXmax();
1299                                                 Int_t KaonMinTOF = (fDiHadronTPCTOF[1][ptbin]->GetAxis(3))->GetXmin();
1300                                                 Int_t KaonMaxTOF = (fDiHadronTPCTOF[1][ptbin]->GetAxis(3))->GetXmax();
1301                                                 if (KaonFillTPC<KaonMaxTPC&&KaonFillTPC>KaonMinTPC&&KaonFillTOF<KaonMaxTOF&&KaonFillTOF>KaonMinTOF) {
1302                                                         fMixedEventsTPCTOFCut[1]->Fill(DPhi,DEta,pt);
1303                                                 }
1304
1305                                                 Double_t ProtonFillTPC = TPCsignal - expectedTPCsignalProton;
1306                                                 Double_t ProtonFillTOF = TOFsignal - expectedTOFsignalProton;
1307                                                 Int_t ProtonMinTPC = (fDiHadronTPCTOF[2][ptbin]->GetAxis(2))->GetXmin();
1308                                                 Int_t ProtonMaxTPC = (fDiHadronTPCTOF[2][ptbin]->GetAxis(2))->GetXmax();
1309                                                 Int_t ProtonMinTOF = (fDiHadronTPCTOF[2][ptbin]->GetAxis(3))->GetXmin();
1310                                                 Int_t ProtonMaxTOF = (fDiHadronTPCTOF[2][ptbin]->GetAxis(3))->GetXmax();
1311                                                 if (ProtonFillTPC<ProtonMaxTPC&&ProtonFillTPC>ProtonMinTPC&&ProtonFillTOF<ProtonMaxTOF&&ProtonFillTOF>ProtonMinTOF) {
1312                                                         fMixedEventsTPCTOFCut[2]->Fill(DPhi,DEta,pt);
1313                                                 }
1314                                                 */
1315
1316                     }
1317                 }
1318             }
1319         }
1320
1321         // Copy the triggers from the current event into the buffer.
1322         if (fAODVertex->GetZ()<10.) {
1323
1324             if (fVerbose>3) cout<<"AliAnalysisTaskDiHadronPID::UserExec -> Copying "<<triggers->GetEntriesFast()<<" triggers with vertex z = "<<fAODVertex->GetZ()<<" to the buffer."<<endl;
1325
1326             for (Int_t iTrig = 0; iTrig<triggers->GetEntriesFast(); iTrig++) {
1327
1328                 currentTrigger = (AliAODTrack*)(triggers->At(iTrig));
1329                 if (fVerbose>3) cout<<"AliAnalysisTaskDiHadronPID::UserExec -> Trigger pt = "<<currentTrigger->Pt()<<endl;
1330
1331                 fTrigBuffer[fTrigBufferIndex][0] = fAODVertex->GetZ();
1332                 fTrigBuffer[fTrigBufferIndex][1] = currentTrigger->Phi();
1333                 fTrigBuffer[fTrigBufferIndex][2] = currentTrigger->Eta();
1334                 fTrigBuffer[fTrigBufferIndex][3] = currentTrigger->Pt();
1335                 fTrigBufferIndex++;
1336                 if (fTrigBufferSize<fTrigBufferMaxSize) {fTrigBufferSize++;} // 250 triggers should be enough to get 10 times more data in mixed events.
1337                 if (fTrigBufferIndex==fTrigBufferMaxSize) {fTrigBufferIndex=0;}
1338             }
1339         }
1340     }
1341
1342     if (fPrintBufferSize) cout<<"AliAnalysisTaskDiHadronPID::UserExec -> Trigger buffer index: "<<fTrigBufferIndex<<", and size: "<<fTrigBufferSize<<endl;
1343
1344         delete triggers;
1345         delete associateds;
1346
1347         PostData(1,fHistoList);
1348
1349 }
1350
1351 //_____________________________________________________________________________
1352 void AliAnalysisTaskDiHadronPID::Terminate(Option_t *)
1353
1354 {
1355         //
1356         // Terminate.
1357     //
1358
1359 }
1360