PWGHF/vertexingHF/AliAnalysisTaskCombinHF.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2018, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id: $ */
  17
  18 //*************************************************************************
  19 // Class AliAnalysisTaskCombinHF
  20 // AliAnalysisTaskSE to build D meson candidates by combining tracks
  21 //  background is computed LS and track rotations is
  22 // Authors: F. Prino, A. Rossi
  23 /////////////////////////////////////////////////////////////
  24
  25 #include <TList.h>
  26 #include <TH1F.h>
  27 #include <TH2F.h>
  28 #include <TH3F.h>
  29 #include <THnSparse.h>
  30
  31 #include "AliAnalysisManager.h"
  32 #include "AliInputEventHandler.h"
  33 #include "AliPIDResponse.h"
  34 #include "AliAODHandler.h"
  35 #include "AliAODEvent.h"
  36 #include "AliAODMCParticle.h"
  37 #include "AliAODMCHeader.h"
  38 #include "AliAODVertex.h"
  39 #include "AliAODTrack.h"
  40 #include "AliVertexingHFUtils.h"
  41 #include "AliAnalysisTaskCombinHF.h"
  42
  43 ClassImp(AliAnalysisTaskCombinHF)
  44
  45
  46 //________________________________________________________________________
  47 AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF():
  48   AliAnalysisTaskSE(),
  49   fOutput(0x0),
  50   fHistNEvents(0x0),
  51   fHistTrackStatus(0x0),
  52   fHistCheckOrigin(0x0),
  53   fHistCheckOriginSel(0x0),
  54   fHistCheckDecChan(0x0),
  55   fHistCheckDecChanAcc(0x0),
  56   fPtVsYVsMultGen(0x0),
  57   fPtVsYVsMultGenLargeAcc(0x0),
  58   fPtVsYVsMultGenLimAcc(0x0),
  59   fPtVsYVsMultGenAcc(0x0),
  60   fPtVsYVsMultReco(0x0),
  61   fMassVsPtVsY(0x0),
  62   fMassVsPtVsYRot(0x0),
  63   fMassVsPtVsYLSpp(0x0),
  64   fMassVsPtVsYLSmm(0x0),
  65   fMassVsPtVsYSig(0x0),
  66   fMassVsPtVsYRefl(0x0),
  67   fMassVsPtVsYBkg(0x0),
  68   fNSelected(0x0),
  69   fNormRotated(0x0),
  70   fDeltaMass(0x0),
  71   fDeltaMassFullAnalysis(0x0),
  72   fMassVsPtVsYME(0x0),
  73   fMassVsPtVsYMELSpp(0x0),
  74   fMassVsPtVsYMELSmm(0x0),
  75   fEventsPerPool(0x0),
  76   fMixingsPerPool(0x0),
  77   fFilterMask(BIT(4)),
  78   fTrackCutsAll(0x0),
  79   fTrackCutsPion(0x0),
  80   fTrackCutsKaon(0x0),
  81   fPidHF(new AliAODPidHF()),
  82   fAnalysisCuts(0x0),
  83   fMinMass(1.720),
  84   fMaxMass(2.150),
  85   fMaxPt(10.),
  86   fPtBinWidth(0.5),
  87   fEtaAccCut(0.9),
  88   fPtAccCut(0.1),
  89   fNRotations(9),
  90   fMinAngleForRot(5*TMath::Pi()/6),
  91   fMaxAngleForRot(7*TMath::Pi()/6),
  92   fMinAngleForRot3(2*TMath::Pi()/6),
  93   fMaxAngleForRot3(4*TMath::Pi()/6),
  94   fCounter(0x0),
  95   fMeson(kDzero),
  96   fReadMC(kFALSE),
  97   fPromptFeeddown(kPrompt),
  98   fGoUpToQuark(kTRUE),
  99   fFullAnalysis(0),
 100   fPIDstrategy(knSigma),
 101   fmaxPforIDPion(0.8),
 102   fmaxPforIDKaon(2.),
 103   fKeepNegID(kFALSE),
 104   fPIDselCaseZero(0),
 105   fBayesThresKaon(0.4),
 106   fBayesThresPion(0.4),
 107   fDoEventMixing(1),
 108   fNumberOfEventsForMixing(20),
 109   fMaxzVertDistForMix(5.),
 110   fMaxMultDiffForMix(5.),
 111   fNzVertPools(1),
 112   fNzVertPoolsLimSize(2),
 113   fzVertPoolLims(0x0),
 114   fNMultPools(1),
 115   fNMultPoolsLimSize(2),
 116   fMultPoolLims(0x0),
 117   fNOfPools(1),
 118   fEventBuffer(0x0),
 119   fEventInfo(new TObjString("")),
 120   fVtxZ(0),
 121   fMultiplicity(0),
 122   fMinMultiplicity(-0.5),
 123   fMaxMultiplicity(199.5),
 124   fKaonTracks(0x0),
 125   fPionTracks(0x0)
 126 {
 127   // default constructor
 128 }
 129
 130 //________________________________________________________________________
 131 AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF(Int_t meson, AliRDHFCuts* analysiscuts):
 132   AliAnalysisTaskSE("DmesonCombin"),
 133   fOutput(0x0),
 134   fHistNEvents(0x0),
 135   fHistTrackStatus(0x0),
 136   fHistCheckOrigin(0x0),
 137   fHistCheckOriginSel(0x0),
 138   fHistCheckDecChan(0x0),
 139   fHistCheckDecChanAcc(0x0),
 140   fPtVsYVsMultGen(0x0),
 141   fPtVsYVsMultGenLargeAcc(0x0),
 142   fPtVsYVsMultGenLimAcc(0x0),
 143   fPtVsYVsMultGenAcc(0x0),
 144   fPtVsYVsMultReco(0x0),
 145   fMassVsPtVsY(0x0),
 146   fMassVsPtVsYRot(0x0),
 147   fMassVsPtVsYLSpp(0x0),
 148   fMassVsPtVsYLSmm(0x0),
 149   fMassVsPtVsYSig(0x0),
 150   fMassVsPtVsYRefl(0x0),
 151   fMassVsPtVsYBkg(0x0),
 152   fNSelected(0x0),
 153   fNormRotated(0x0),
 154   fDeltaMass(0x0),
 155   fDeltaMassFullAnalysis(0x0),
 156   fMassVsPtVsYME(0x0),
 157   fMassVsPtVsYMELSpp(0x0),
 158   fMassVsPtVsYMELSmm(0x0),
 159   fEventsPerPool(0x0),
 160   fMixingsPerPool(0x0),
 161   fFilterMask(BIT(4)),
 162   fTrackCutsAll(0x0),
 163   fTrackCutsPion(0x0),
 164   fTrackCutsKaon(0x0),
 165   fPidHF(new AliAODPidHF()),
 166   fAnalysisCuts(analysiscuts),
 167   fMinMass(1.720),
 168   fMaxMass(2.150),
 169   fMaxPt(10.),
 170   fPtBinWidth(0.5),
 171   fEtaAccCut(0.9),
 172   fPtAccCut(0.1),
 173   fNRotations(9),
 174   fMinAngleForRot(5*TMath::Pi()/6),
 175   fMaxAngleForRot(7*TMath::Pi()/6),
 176   fMinAngleForRot3(2*TMath::Pi()/6),
 177   fMaxAngleForRot3(4*TMath::Pi()/6),
 178   fCounter(0x0),
 179   fMeson(meson),
 180   fReadMC(kFALSE),
 181   fPromptFeeddown(1),
 182   fGoUpToQuark(kTRUE),
 183   fFullAnalysis(0),
 184   fPIDstrategy(knSigma),
 185   fmaxPforIDPion(0.8),
 186   fmaxPforIDKaon(2.),
 187   fKeepNegID(kFALSE),
 188   fPIDselCaseZero(0),
 189   fBayesThresKaon(0.4),
 190   fBayesThresPion(0.4),
 191   fDoEventMixing(1),
 192   fNumberOfEventsForMixing(20),
 193   fMaxzVertDistForMix(5.),
 194   fMaxMultDiffForMix(5.),
 195   fNzVertPools(1),
 196   fNzVertPoolsLimSize(2),
 197   fzVertPoolLims(0x0),
 198   fNMultPools(1),
 199   fNMultPoolsLimSize(2),
 200   fMultPoolLims(0x0),
 201   fNOfPools(1),
 202   fEventBuffer(0x0),
 203   fEventInfo(new TObjString("")),
 204   fVtxZ(0),
 205   fMultiplicity(0),
 206   fMinMultiplicity(-0.5),
 207   fMaxMultiplicity(199.5),
 208   fKaonTracks(0x0),
 209   fPionTracks(0x0)
 210 {
 211   // standard constructor
 212
 213   DefineOutput(1,TList::Class());  //My private output
 214   DefineOutput(2,AliNormalizationCounter::Class());
 215 }
 216
 217 //________________________________________________________________________
 218 AliAnalysisTaskCombinHF::~AliAnalysisTaskCombinHF()
 219 {
 220   //
 221   // Destructor
 222   //
 223   if(fOutput && !fOutput->IsOwner()){
 224     delete fHistNEvents;
 225     delete fHistTrackStatus;
 226     delete fHistCheckOrigin;
 227     delete fHistCheckOriginSel;
 228     delete fHistCheckDecChan;
 229     delete fHistCheckDecChanAcc;
 230     delete fPtVsYVsMultGen;
 231     delete fPtVsYVsMultGenLargeAcc;
 232     delete fPtVsYVsMultGenLimAcc;
 233     delete fPtVsYVsMultGenAcc;
 234     delete fPtVsYVsMultReco;
 235     delete fMassVsPtVsY;
 236     delete fMassVsPtVsYLSpp;
 237     delete fMassVsPtVsYLSmm;
 238     delete fMassVsPtVsYRot;
 239     delete fMassVsPtVsYSig;
 240     delete fMassVsPtVsYRefl;
 241     delete fMassVsPtVsYBkg;
 242     delete fNSelected;
 243     delete fNormRotated;
 244     delete fDeltaMass;
 245     delete fDeltaMassFullAnalysis;
 246     delete fMassVsPtVsYME;
 247     delete fMassVsPtVsYMELSpp;
 248     delete fMassVsPtVsYMELSmm;
 249   }
 250
 251   delete fOutput;
 252   delete fCounter;
 253   delete fTrackCutsAll;
 254   delete fTrackCutsPion;
 255   delete fTrackCutsKaon;
 256   delete fPidHF;
 257   delete fAnalysisCuts;
 258   if(fKaonTracks) fKaonTracks->Delete();
 259   if(fPionTracks) fPionTracks->Delete();
 260   delete fKaonTracks;
 261   delete fPionTracks;
 262
 263   if(fEventBuffer){
 264     for(Int_t i=0; i<fNOfPools; i++) delete fEventBuffer[i];
 265     delete fEventBuffer;
 266   }
 267   delete fEventInfo;
 268   delete [] fzVertPoolLims;
 269   delete [] fMultPoolLims;
 270 }
 271
 272 //________________________________________________________________________
 273 void AliAnalysisTaskCombinHF::ConfigureZVertPools(Int_t nPools, Double_t*  zVertLimits)
 274 {
 275   // sets the pools for event mizing in zvertex
 276   if(fzVertPoolLims) delete [] fzVertPoolLims;
 277   fNzVertPools=nPools;
 278   fNzVertPoolsLimSize=nPools+1;
 279   fzVertPoolLims = new Double_t[fNzVertPoolsLimSize];
 280   for(Int_t ib=0; ib<fNzVertPoolsLimSize; ib++) fzVertPoolLims[ib]=zVertLimits[ib];
 281   return;
 282 }
 283 //________________________________________________________________________
 284 void AliAnalysisTaskCombinHF::ConfigureMultiplicityPools(Int_t nPools, Double_t*  multLimits)
 285 {
 286   // sets the pools for event mizing in zvertex
 287   if(fMultPoolLims) delete [] fMultPoolLims;
 288   fNMultPools=nPools;
 289   fNMultPoolsLimSize=nPools+1;
 290   fMultPoolLims = new Double_t[fNMultPoolsLimSize];
 291   for(Int_t ib=0; ib<nPools+1; ib++) fMultPoolLims[ib]=multLimits[ib];
 292   return;
 293 }
 294 //________________________________________________________________________
 295 void AliAnalysisTaskCombinHF::UserCreateOutputObjects()
 296 {
 297   // Create the output container
 298   //
 299   if(fDebug > 1) printf("AnalysisTaskCombinHF::UserCreateOutputObjects() \n");
 300
 301   fOutput = new TList();
 302   fOutput->SetOwner();
 303   fOutput->SetName("OutputHistos");
 304
 305   fHistNEvents = new TH1F("hNEvents", "number of events ",9,-0.5,8.5);
 306   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 307   fHistNEvents->GetXaxis()->SetBinLabel(2,"n. passing IsEvSelected");
 308   fHistNEvents->GetXaxis()->SetBinLabel(3,"n. rejected due to trigger");
 309   fHistNEvents->GetXaxis()->SetBinLabel(4,"n. rejected due to phys sel");
 310   fHistNEvents->GetXaxis()->SetBinLabel(5,"n. rejected due to not reco vertex");
 311   fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected for contr vertex");
 312   fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected for vertex out of accept");
 313   fHistNEvents->GetXaxis()->SetBinLabel(8,"n. rejected for pileup events");
 314   fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of out centrality events");
 315
 316   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
 317   fHistNEvents->Sumw2();
 318   fHistNEvents->SetMinimum(0);
 319   fOutput->Add(fHistNEvents);
 320
 321   fHistTrackStatus  = new TH1F("hTrackStatus", "",8,-0.5,7.5);
 322   fHistTrackStatus->GetXaxis()->SetBinLabel(1,"Not OK");
 323   fHistTrackStatus->GetXaxis()->SetBinLabel(2,"Track OK");
 324   fHistTrackStatus->GetXaxis()->SetBinLabel(3,"Kaon, Not OK");
 325   fHistTrackStatus->GetXaxis()->SetBinLabel(4,"Kaon OK");
 326   fHistTrackStatus->GetXaxis()->SetBinLabel(5,"Pion, Not OK");
 327   fHistTrackStatus->GetXaxis()->SetBinLabel(6,"Pion OK");
 328   fHistTrackStatus->GetXaxis()->SetBinLabel(7,"Kaon||Pion, Not OK");
 329   fHistTrackStatus->GetXaxis()->SetBinLabel(8,"Kaon||Pion OK");
 330   fHistTrackStatus->GetXaxis()->SetNdivisions(1,kFALSE);
 331   fHistTrackStatus->Sumw2();
 332   fHistTrackStatus->SetMinimum(0);
 333   fOutput->Add(fHistTrackStatus);
 334
 335   Int_t nPtBins = (Int_t)(fMaxPt/fPtBinWidth+0.001);
 336   Double_t maxPt=fPtBinWidth*nPtBins;
 337
 338   if(fReadMC){
 339
 340     fHistCheckOrigin=new TH1F("hCheckOrigin","",7,-1.5,5.5);
 341     fHistCheckOrigin->Sumw2();
 342     fHistCheckOrigin->SetMinimum(0);
 343     fOutput->Add(fHistCheckOrigin);
 344
 345     fHistCheckOriginSel=new TH1F("hCheckOriginSel","",7,-1.5,5.5);
 346     fHistCheckOriginSel->Sumw2();
 347     fHistCheckOriginSel->SetMinimum(0);
 348     fOutput->Add(fHistCheckOriginSel);
 349
 350     fHistCheckDecChan=new TH1F("hCheckDecChan","",7,-2.5,4.5);
 351     fHistCheckDecChan->Sumw2();
 352     fHistCheckDecChan->SetMinimum(0);
 353     fOutput->Add(fHistCheckDecChan);
 354
 355     fHistCheckDecChanAcc=new TH1F("hCheckDecChanAcc","",7,-2.5,4.5);
 356     fHistCheckDecChanAcc->Sumw2();
 357     fHistCheckDecChanAcc->SetMinimum(0);
 358     fOutput->Add(fHistCheckDecChanAcc);
 359
 360     fPtVsYVsMultGen= new TH3F("hPtVsYVsMultGen","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
 361     fPtVsYVsMultGen->Sumw2();
 362     fPtVsYVsMultGen->SetMinimum(0);
 363     fOutput->Add(fPtVsYVsMultGen);
 364
 365     fPtVsYVsMultGenLargeAcc= new TH3F("hPtVsYVsMultGenLargeAcc","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
 366     fPtVsYVsMultGenLargeAcc->Sumw2();
 367     fPtVsYVsMultGenLargeAcc->SetMinimum(0);
 368     fOutput->Add(fPtVsYVsMultGenLargeAcc);
 369
 370     fPtVsYVsMultGenLimAcc= new TH3F("hPtVsYVsMultGenLimAcc","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
 371     fPtVsYVsMultGenLimAcc->Sumw2();
 372     fPtVsYVsMultGenLimAcc->SetMinimum(0);
 373     fOutput->Add(fPtVsYVsMultGenLimAcc);
 374
 375     fPtVsYVsMultGenAcc= new TH3F("hPtVsYVsMultGenAcc","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
 376     fPtVsYVsMultGenAcc->Sumw2();
 377     fPtVsYVsMultGenAcc->SetMinimum(0);
 378     fOutput->Add(fPtVsYVsMultGenAcc);
 379
 380     fPtVsYVsMultReco= new TH3F("hPtVsYVsMultReco","",nPtBins,0.,maxPt,20,-1.,1.,200,fMinMultiplicity,fMaxMultiplicity);
 381     fPtVsYVsMultReco->Sumw2();
 382     fPtVsYVsMultReco->SetMinimum(0);
 383     fOutput->Add(fPtVsYVsMultReco);
 384   }
 385
 386
 387   Int_t nMassBins=static_cast<Int_t>(fMaxMass*1000.-fMinMass*1000.);
 388   Double_t maxm=fMinMass+nMassBins*0.001;
 389   fMassVsPtVsY=new TH3F("hMassVsPtVsY","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
 390   fMassVsPtVsY->Sumw2();
 391   fMassVsPtVsY->SetMinimum(0);
 392   fOutput->Add(fMassVsPtVsY);
 393
 394   fMassVsPtVsYRot=new TH3F("hMassVsPtVsYRot","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
 395   fMassVsPtVsYRot->Sumw2();
 396   fMassVsPtVsYRot->SetMinimum(0);
 397   fOutput->Add(fMassVsPtVsYRot);
 398
 399   if(fMeson==kDzero){
 400     fMassVsPtVsYLSpp=new TH3F("hMassVsPtVsYLSpp","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
 401     fMassVsPtVsYLSpp->Sumw2();
 402     fMassVsPtVsYLSpp->SetMinimum(0);
 403     fOutput->Add(fMassVsPtVsYLSpp);
 404     fMassVsPtVsYLSmm=new TH3F("hMassVsPtVsYLSmm","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
 405     fMassVsPtVsYLSmm->Sumw2();
 406     fMassVsPtVsYLSmm->SetMinimum(0);
 407     fOutput->Add(fMassVsPtVsYLSmm);
 408   }
 409
 410   fMassVsPtVsYSig=new TH3F("hMassVsPtVsYSig","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
 411   fMassVsPtVsYSig->Sumw2();
 412   fMassVsPtVsYSig->SetMinimum(0);
 413   fOutput->Add(fMassVsPtVsYSig);
 414
 415   fMassVsPtVsYRefl=new TH3F("hMassVsPtVsYRefl","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
 416   fMassVsPtVsYRefl->Sumw2();
 417   fMassVsPtVsYRefl->SetMinimum(0);
 418   fOutput->Add(fMassVsPtVsYRefl);
 419
 420   fMassVsPtVsYBkg=new TH3F("hMassVsPtVsYBkg","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
 421   fMassVsPtVsYBkg->Sumw2();
 422   fMassVsPtVsYBkg->SetMinimum(0);
 423   fOutput->Add(fMassVsPtVsYBkg);
 424
 425   fNSelected=new TH1F("hNSelected","",100,-0.5,99.5);
 426   fNSelected->Sumw2();
 427   fNSelected->SetMinimum(0);
 428   fOutput->Add(fNSelected);
 429
 430   fNormRotated=new TH1F("hNormRotated","",11,-0.5,10.5);
 431   fNormRotated->Sumw2();
 432   fNormRotated->SetMinimum(0);
 433   fOutput->Add(fNormRotated);
 434
 435   fDeltaMass=new TH1F("hDeltaMass","",100,-0.4,0.4);
 436   fDeltaMass->Sumw2();
 437   fDeltaMass->SetMinimum(0);
 438   fOutput->Add(fDeltaMass);
 439
 440   Int_t binSparseDMassRot[5]={nMassBins,100,24,40,20};
 441   Double_t edgeLowSparseDMassRot[5]={fMinMass,-0.4,0.,-4.,0};
 442   Double_t edgeHighSparseDMassRot[5]={maxm,0.4,12.,4.,3.14};
 443   fDeltaMassFullAnalysis=new THnSparseF("fDeltaMassFullAnalysis","fDeltaMassFullAnalysis;inv mass (GeV/c);#Delta inv mass (GeV/c) ; p_{T}^{D} (GeV/c); #Delta p_{T} (GeV/c); daughter angle (2prongs) (rad);",5,binSparseDMassRot,edgeLowSparseDMassRot,edgeHighSparseDMassRot);
 444   fOutput->Add(fDeltaMassFullAnalysis);
 445
 446   fMassVsPtVsYME=new TH3F("hMassVsPtVsYME","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
 447   fMassVsPtVsYME->Sumw2();
 448   fMassVsPtVsYME->SetMinimum(0);
 449   fOutput->Add(fMassVsPtVsYME);
 450
 451   fMassVsPtVsYMELSpp=new TH3F("hMassVsPtVsYMELSpp","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
 452   fMassVsPtVsYMELSpp->Sumw2();
 453   fMassVsPtVsYMELSpp->SetMinimum(0);
 454   fOutput->Add(fMassVsPtVsYMELSpp);
 455
 456   fMassVsPtVsYMELSmm=new TH3F("hMassVsPtVsYMELSmm","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
 457   fMassVsPtVsYMELSmm->Sumw2();
 458   fMassVsPtVsYMELSmm->SetMinimum(0);
 459   fOutput->Add(fMassVsPtVsYMELSmm);
 460
 461   fNOfPools=fNzVertPools*fNMultPools;
 462   if(!fzVertPoolLims || !fMultPoolLims) fNOfPools=1;
 463   if(fDoEventMixing==2) fNOfPools=1;
 464   if(fNOfPools>1 && fzVertPoolLims && fMultPoolLims){
 465     fEventsPerPool=new TH2F("hEventsPerPool","hEventsPerPool",fNzVertPools,fzVertPoolLims,fNMultPools,fMultPoolLims);
 466     fMixingsPerPool=new TH2F("hMixingsPerPool","hMixingsPerPool",fNzVertPools,fzVertPoolLims,fNMultPools,fMultPoolLims);
 467   }else{
 468     fEventsPerPool=new TH2F("hEventsPerPool","hEventsPerPool",1,-10.,10.,1,-0.5,2000.5);
 469     fMixingsPerPool=new TH2F("hMixingsPerPool","hMixingsPerPool",1,-10.,10.,1,-0.5,2000.5);
 470   }
 471   fEventsPerPool->Sumw2();
 472   fEventsPerPool->SetMinimum(0);
 473   fOutput->Add(fEventsPerPool);
 474   fMixingsPerPool->Sumw2();
 475   fMixingsPerPool->SetMinimum(0);
 476   fOutput->Add(fMixingsPerPool);
 477
 478   //Counter for Normalization
 479   fCounter = new AliNormalizationCounter("NormalizationCounter");
 480   fCounter->Init();
 481
 482   fKaonTracks = new TObjArray();
 483   fPionTracks=new TObjArray();
 484   fKaonTracks->SetOwner();
 485   fPionTracks->SetOwner();
 486
 487   fEventBuffer = new TTree*[fNOfPools];
 488   for(Int_t i=0; i<fNOfPools; i++){
 489     fEventBuffer[i]=new TTree(Form("EventBuffer_%d",i), "Temporary buffer for event mixing");
 490     fEventBuffer[i]->Branch("zVertex", &fVtxZ);
 491     fEventBuffer[i]->Branch("multiplicity", &fMultiplicity);
 492     fEventBuffer[i]->Branch("eventInfo", "TObjString",&fEventInfo);
 493     fEventBuffer[i]->Branch("karray", "TObjArray", &fKaonTracks);
 494     fEventBuffer[i]->Branch("parray", "TObjArray", &fPionTracks);
 495   }
 496
 497   PostData(1,fOutput);
 498   PostData(2,fCounter);
 499 }
 500
 501 //________________________________________________________________________
 502 void AliAnalysisTaskCombinHF::UserExec(Option_t */*option*/){
 503   //Build the 3-track combinatorics (+-+ and -+-) for D+->Kpipi decays
 504
 505   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 506   if(!aod && AODEvent() && IsStandardAOD()) {
 507     // In case there is an AOD handler writing a standard AOD, use the AOD
 508     // event in memory rather than the input (ESD) event.
 509     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 510   }
 511   if(!aod){
 512     printf("AliAnalysisTaskCombinHF::UserExec: AOD not found!\n");
 513     return;
 514   }
 515
 516   // fix for temporary bug in ESDfilter
 517   // the AODs with null vertex pointer didn't pass the PhysSel
 518   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
 519   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 520   AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
 521   AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
 522   fPidHF->SetPidResponse(pidResp);
 523
 524
 525   fHistNEvents->Fill(0); // count event
 526   // Post the data already here
 527   PostData(1,fOutput);
 528
 529   fCounter->StoreEvent(aod,fAnalysisCuts,fReadMC);
 530
 531   Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
 532   if(fAnalysisCuts->IsEventRejectedDueToTrigger() || fAnalysisCuts->IsEventRejectedDuePhysicsSelection()){
 533     if(fAnalysisCuts->IsEventRejectedDueToTrigger()) fHistNEvents->Fill(2);
 534     if(fAnalysisCuts->IsEventRejectedDuePhysicsSelection()) fHistNEvents->Fill(3);
 535   }else{
 536     if(fAnalysisCuts->IsEventRejectedDueToCentrality()){
 537       fHistNEvents->Fill(8);
 538     }else{
 539       if(fAnalysisCuts->IsEventRejectedDueToNotRecoVertex() || fAnalysisCuts->IsEventRejectedDueToVertexContributors()){
 540         if(fAnalysisCuts->IsEventRejectedDueToNotRecoVertex())fHistNEvents->Fill(4);
 541         if(fAnalysisCuts->IsEventRejectedDueToVertexContributors())fHistNEvents->Fill(5);
 542       }else{
 543         if(fAnalysisCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion())fHistNEvents->Fill(6);
 544         else if(fAnalysisCuts->IsEventRejectedDueToPileup())fHistNEvents->Fill(7);
 545       }
 546     }
 547   }
 548
 549
 550   if(!isEvSel)return;
 551
 552   fHistNEvents->Fill(1);
 553
 554   TClonesArray *arrayMC=0;
 555   AliAODMCHeader *mcHeader=0;
 556   if(fReadMC){
 557     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 558     if(!arrayMC) {
 559       printf("AliAnalysisTaskCombinHF::UserExec: MC particles branch not found!\n");
 560       return;
 561     }
 562
 563     // load MC header
 564     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 565     if(!mcHeader) {
 566       printf("AliAnalysisTaskCombinHF::UserExec: MC header branch not found!\n");
 567       return;
 568     }
 569     FillGenHistos(arrayMC);
 570   }
 571
 572   Int_t ntracks=aod->GetNumberOfTracks();
 573   fVtxZ = aod->GetPrimaryVertex()->GetZ();
 574   fMultiplicity = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
 575
 576   // select and flag tracks
 577   UChar_t* status = new UChar_t[ntracks];
 578   for(Int_t iTr=0; iTr<ntracks; iTr++){
 579     status[iTr]=0;
 580     AliAODTrack* track=dynamic_cast<AliAODTrack*>(aod->GetTrack(iTr));
 581     if(!track) AliFatal("Not a standard AOD");
 582     if(IsTrackSelected(track)) status[iTr]+=1;
 583
 584     // PID
 585     if (fPIDstrategy == knSigma) {
 586       // nsigma PID
 587       if(IsKaon(track)) status[iTr]+=2;
 588       if(IsPion(track)) status[iTr]+=4;
 589     }
 590     else if (fPIDstrategy == kBayesianMaxProb || fPIDstrategy == kBayesianThres) {
 591       // Bayesian PID
 592       Double_t *weights = new Double_t[AliPID::kSPECIES];
 593       fPidHF->GetPidCombined()->ComputeProbabilities(track, fPidHF->GetPidResponse(), weights);
 594       if (fPIDstrategy == kBayesianMaxProb) {
 595         if (TMath::MaxElement(AliPID::kSPECIES, weights) == weights[AliPID::kKaon]) status[iTr] += 2;
 596         if (TMath::MaxElement(AliPID::kSPECIES, weights) == weights[AliPID::kPion]) status[iTr] += 4;
 597       }
 598       if (fPIDstrategy == kBayesianThres) {
 599         if (weights[AliPID::kKaon] > fBayesThresKaon) status[iTr] += 2;
 600         if (weights[AliPID::kPion] > fBayesThresPion) status[iTr] += 4;
 601       }
 602       delete[] weights;
 603     }
 604
 605     fHistTrackStatus->Fill(status[iTr]);
 606   }
 607
 608   // build the combinatorics
 609   Int_t nSelected=0;
 610   Int_t nFiltered=0;
 611   Double_t dummypos[3]={0.,0.,0.};
 612   AliAODVertex* v2=new AliAODVertex(dummypos,999.,-1,2);
 613   AliAODVertex* v3=new AliAODVertex(dummypos,999.,-1,3);
 614   // dummy values of track impact parameter, needed to build an AliAODRecoDecay object
 615   Double_t d02[2]={0.,0.};
 616   Double_t d03[3]={0.,0.,0.};
 617   AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
 618   AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
 619   UInt_t pdg0[2]={321,211};
 620   UInt_t pdgp[3]={321,211,211};
 621   //  UInt_t pdgs[3]={321,321,211};
 622   Double_t tmpp[3];
 623   Double_t px[3],py[3],pz[3];
 624   Int_t dgLabels[3];
 625   fKaonTracks->Delete();
 626   fPionTracks->Delete();
 627
 628   for(Int_t iTr1=0; iTr1<ntracks; iTr1++){
 629     AliAODTrack* trK=dynamic_cast<AliAODTrack*>(aod->GetTrack(iTr1));
 630     if(!trK) AliFatal("Not a standard AOD");
 631     if((status[iTr1] & 1)==0) continue;
 632     if(fDoEventMixing>0){
 633       if(status[iTr1] & 2) fKaonTracks->AddLast(new TLorentzVector(trK->Px(),trK->Py(),trK->Pz(),trK->Charge()));
 634       if(status[iTr1] & 4) fPionTracks->AddLast(new TLorentzVector(trK->Px(),trK->Py(),trK->Pz(),trK->Charge()));
 635     }
 636     if((status[iTr1] & 2)==0) continue;
 637     Int_t chargeK=trK->Charge();
 638     trK->GetPxPyPz(tmpp);
 639     px[0] = tmpp[0];
 640     py[0] = tmpp[1];
 641     pz[0] = tmpp[2];
 642     dgLabels[0]=trK->GetLabel();
 643     for(Int_t iTr2=0; iTr2<ntracks; iTr2++){
 644       if((status[iTr2] & 1)==0) continue;
 645       if((status[iTr2] & 4)==0) continue;
 646       if(iTr1==iTr2) continue;
 647       AliAODTrack* trPi1=dynamic_cast<AliAODTrack*>(aod->GetTrack(iTr2));
 648       if(!trPi1) AliFatal("Not a standard AOD");
 649       Int_t chargePi1=trPi1->Charge();
 650       trPi1->GetPxPyPz(tmpp);
 651       px[1] = tmpp[0];
 652       py[1] = tmpp[1];
 653       pz[1] = tmpp[2];
 654       dgLabels[1]=trPi1->GetLabel();
 655       if(chargePi1==chargeK){
 656         if(fMeson==kDzero) FillLSHistos(421,2,tmpRD2,px,py,pz,pdg0,chargePi1);
 657         continue;
 658       }
 659       if(fMeson==kDzero){
 660         nFiltered++;
 661         v2->AddDaughter(trK);
 662         v2->AddDaughter(trPi1);
 663         tmpRD2->SetSecondaryVtx(v2);
 664         Bool_t ok=FillHistos(421,2,tmpRD2,px,py,pz,pdg0,arrayMC,dgLabels);
 665         v2->RemoveDaughters();
 666         if(ok) nSelected++;
 667       }else{
 668         for(Int_t iTr3=iTr2+1; iTr3<ntracks; iTr3++){
 669           if((status[iTr3] & 1)==0) continue;
 670           if((status[iTr3] & 4)==0) continue;
 671           if(iTr1==iTr3) continue;
 672           AliAODTrack* trPi2=dynamic_cast<AliAODTrack*>(aod->GetTrack(iTr3));
 673           if(!trPi2) AliFatal("Not a standard AOD");
 674           Int_t chargePi2=trPi2->Charge();
 675           if(chargePi2==chargeK) continue;
 676           nFiltered++;
 677           trPi2->GetPxPyPz(tmpp);
 678           px[2] = tmpp[0];
 679           py[2] = tmpp[1];
 680           pz[2] = tmpp[2];
 681           dgLabels[2]=trPi2->GetLabel();
 682           v3->AddDaughter(trK);
 683           v3->AddDaughter(trPi1);
 684           v3->AddDaughter(trPi2);
 685           tmpRD3->SetSecondaryVtx(v3);
 686           Bool_t ok=FillHistos(411,3,tmpRD3,px,py,pz,pdgp,arrayMC,dgLabels);
 687           v3->RemoveDaughters();
 688           if(ok) nSelected++;
 689         }
 690       }
 691     }
 692   }
 693
 694   delete [] status;
 695   delete v2;
 696   delete v3;
 697   delete tmpRD2;
 698   delete tmpRD3;
 699
 700   fNSelected->Fill(nSelected);
 701
 702   fCounter->StoreCandidates(aod,nFiltered,kTRUE);
 703   fCounter->StoreCandidates(aod,nSelected,kFALSE);
 704   fEventInfo->SetString(Form("Ev%d_esd%d_Pi%d_K%d",mgr->GetNcalls(),((AliAODHeader*)aod->GetHeader())->GetEventNumberESDFile(),fPionTracks->GetEntries(),fKaonTracks->GetEntries()));
 705   if(fDoEventMixing==1){
 706     Int_t ind=GetPoolIndex(fVtxZ,fMultiplicity);
 707     if(ind>=0 && ind<fNOfPools){
 708       fEventsPerPool->Fill(fVtxZ,fMultiplicity);
 709       fEventBuffer[ind]->Fill();
 710       if(fEventBuffer[ind]->GetEntries() >= fNumberOfEventsForMixing){
 711         fMixingsPerPool->Fill(fVtxZ,fMultiplicity);
 712           DoMixingWithPools(ind);
 713           ResetPool(ind);
 714       }
 715     }
 716   }else if(fDoEventMixing==2){ // mix with cuts, no pools
 717       fEventBuffer[0]->Fill();
 718   }
 719   PostData(1,fOutput);
 720   PostData(2,fCounter);
 721
 722   return;
 723 }
 724
 725 //________________________________________________________________________
 726 void AliAnalysisTaskCombinHF::FillLSHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, Int_t charge){
 727   // Fill histos for LS candidates
 728
 729   tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
 730   Double_t pt = tmpRD->Pt();
 731   Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
 732   if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
 733     Double_t rapid = tmpRD->Y(pdgD);
 734     if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
 735       if(charge>0) fMassVsPtVsYLSpp->Fill(TMath::Sqrt(minv2),pt,rapid);
 736       else fMassVsPtVsYLSmm->Fill(TMath::Sqrt(minv2),pt,rapid);
 737     }
 738   }
 739   return;
 740 }
 741
 742 //________________________________________________________________________
 743 void AliAnalysisTaskCombinHF::FillGenHistos(TClonesArray* arrayMC){
 744   // Fill histos with generated quantities
 745   Int_t totPart=arrayMC->GetEntriesFast();
 746   Int_t thePDG=411;
 747   Int_t nProng=3;
 748   if(fMeson==kDzero){
 749     thePDG=421;
 750     nProng=2;
 751   }
 752   for(Int_t ip=0; ip<totPart; ip++){
 753     AliAODMCParticle *part = (AliAODMCParticle*)arrayMC->At(ip);
 754     if(TMath::Abs(part->GetPdgCode())==thePDG){
 755       Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,part,fGoUpToQuark);
 756       fHistCheckOrigin->Fill(orig);
 757       if(fPromptFeeddown==kFeeddown && orig!=5) continue;
 758       else if(fPromptFeeddown==kPrompt && orig!=4) continue;
 759       else if(fPromptFeeddown==kBoth && orig<4) continue;
 760       fHistCheckOriginSel->Fill(orig);
 761       Int_t deca=0;
 762       Bool_t isGoodDecay=kFALSE;
 763       Int_t labDau[4]={-1,-1,-1,-1};
 764       if(fMeson==kDzero){
 765         deca=AliVertexingHFUtils::CheckD0Decay(arrayMC,part,labDau);
 766         if(part->GetNDaughters()!=2) continue;
 767         if(deca==1) isGoodDecay=kTRUE;
 768       }else if(fMeson==kDplus){
 769         deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,part,labDau);
 770         if(deca>0) isGoodDecay=kTRUE;
 771       }
 772       fHistCheckDecChan->Fill(deca);
 773       if(labDau[0]==-1){
 774         //      printf(Form("Meson %d Label of daughters not filled correctly -- %d\n",fMeson,isGoodDecay));
 775         continue; //protection against unfilled array of labels
 776       }
 777       Bool_t isInAcc=CheckAcceptance(arrayMC,nProng,labDau);
 778       if(isInAcc) fHistCheckDecChanAcc->Fill(deca);
 779       if(isGoodDecay){
 780         Double_t ptgen=part->Pt();
 781         Double_t ygen=part->Y();
 782         if(fAnalysisCuts->IsInFiducialAcceptance(ptgen,ygen)){
 783           fPtVsYVsMultGen->Fill(ptgen,ygen,fMultiplicity);
 784           if(TMath::Abs(ygen)<0.5) fPtVsYVsMultGenLimAcc->Fill(ptgen,ygen,fMultiplicity);
 785           if(isInAcc) fPtVsYVsMultGenAcc->Fill(ptgen,ygen,fMultiplicity);
 786         }
 787         if(TMath::Abs(ygen)<0.9) fPtVsYVsMultGenLargeAcc->Fill(ptgen,ygen,fMultiplicity);
 788       }
 789     }
 790   }
 791 }
 792
 793 //________________________________________________________________________
 794 Bool_t AliAnalysisTaskCombinHF::FillHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, TClonesArray *arrayMC, Int_t* dgLabels){
 795   // Fill histos for candidates with proper charge sign
 796
 797   Bool_t accept=kFALSE;
 798
 799   tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
 800   Double_t pt = tmpRD->Pt();
 801   Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
 802   Double_t mass=TMath::Sqrt(minv2);
 803
 804   if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
 805     Double_t rapid = tmpRD->Y(pdgD);
 806     if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
 807       fMassVsPtVsY->Fill(mass,pt,rapid);
 808       accept=kTRUE;
 809       if(fReadMC){
 810         Int_t signPdg[3]={0,0,0};
 811         for(Int_t iii=0; iii<nProngs; iii++) signPdg[iii]=pdgdau[iii];
 812         Int_t labD = tmpRD->MatchToMC(pdgD,arrayMC,nProngs,signPdg);
 813         if(labD>=0){
 814           AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(TMath::Abs(dgLabels[0])));
 815           if(part){
 816             Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,part,fGoUpToQuark);
 817             if((fPromptFeeddown==kFeeddown && orig==5)|| (fPromptFeeddown==kPrompt && orig==4) || (fPromptFeeddown==kBoth && orig>=4)) {
 818
 819               Int_t pdgCode = TMath::Abs( part->GetPdgCode() );
 820               if(pdgCode==321){
 821                 fMassVsPtVsYSig->Fill(mass,pt,rapid);
 822                 AliAODMCParticle* dmes =  dynamic_cast<AliAODMCParticle*>(arrayMC->At(labD));
 823                 fPtVsYVsMultReco->Fill(dmes->Pt(),dmes->Y(),fMultiplicity);
 824               }else{
 825                 fMassVsPtVsYRefl->Fill(mass,pt,rapid);
 826               }
 827             }
 828           }
 829         }else{
 830           fMassVsPtVsYBkg->Fill(mass,pt,rapid);
 831         }
 832       }
 833     }
 834   }
 835
 836   Int_t nRotated=0;
 837   Double_t massRot=0;// calculated later only if candidate is acceptable
 838   Double_t angleProngXY;
 839   if(TMath::Abs(pdgD)==421)angleProngXY=TMath::ACos((px[0]*px[1]+py[0]*py[1])/TMath::Sqrt((px[0]*px[0]+py[0]*py[0])*(px[1]*px[1]+py[1]*py[1])));
 840   else {
 841     angleProngXY=TMath::ACos(((px[0]+px[1])*px[2]+(py[0]+py[1])*py[2])/TMath::Sqrt(((px[0]+px[1])*(px[0]+px[1])+(py[0]+py[1])*(py[0]+py[1]))*(px[2]*px[2]+py[2]*py[2])));
 842   }
 843   Double_t ptOrig=pt;
 844
 845
 846   Double_t rotStep=(fMaxAngleForRot-fMinAngleForRot)/(fNRotations-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot
 847   Double_t rotStep3=(fMaxAngleForRot3-fMinAngleForRot3)/(fNRotations-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot
 848
 849   for(Int_t irot=0; irot<fNRotations; irot++){
 850     Double_t phirot=fMinAngleForRot+rotStep*irot;
 851     Double_t tmpx=px[0];
 852     Double_t tmpy=py[0];
 853     Double_t tmpx2=px[2];
 854     Double_t tmpy2=py[2];
 855     px[0]=tmpx*TMath::Cos(phirot)-tmpy*TMath::Sin(phirot);
 856     py[0]=tmpx*TMath::Sin(phirot)+tmpy*TMath::Cos(phirot);
 857     if(pdgD==411 || pdgD==431){
 858       Double_t phirot2=fMinAngleForRot3+rotStep3*irot;
 859       px[2]=tmpx*TMath::Cos(phirot2)-tmpy*TMath::Sin(phirot2);
 860       py[2]=tmpx*TMath::Sin(phirot2)+tmpy*TMath::Cos(phirot2);
 861     }
 862     tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
 863     pt = tmpRD->Pt();
 864     minv2 = tmpRD->InvMass2(nProngs,pdgdau);
 865     if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
 866       Double_t rapid = tmpRD->Y(pdgD);
 867       if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
 868         massRot=TMath::Sqrt(minv2);
 869         fMassVsPtVsYRot->Fill(massRot,pt,rapid);
 870         nRotated++;
 871         fDeltaMass->Fill(massRot-mass);
 872         if(fFullAnalysis){
 873           Double_t pointRot[5]={mass,massRot-mass,ptOrig,pt-ptOrig,angleProngXY};
 874           fDeltaMassFullAnalysis->Fill(pointRot);
 875         }
 876       }
 877     }
 878     px[0]=tmpx;
 879     py[0]=tmpy;
 880     if(pdgD==411 || pdgD==431){
 881       px[2]=tmpx2;
 882       py[2]=tmpy2;
 883     }
 884   }
 885   fNormRotated->Fill(nRotated);
 886
 887   return accept;
 888
 889 }
 890 //________________________________________________________________________
 891 void AliAnalysisTaskCombinHF::FillMEHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau){
 892   // Fill histos for candidates in MixedEvents
 893
 894   tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
 895   Double_t pt = tmpRD->Pt();
 896   Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
 897   Double_t mass=TMath::Sqrt(minv2);
 898
 899   if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
 900     Double_t rapid = tmpRD->Y(pdgD);
 901     if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
 902       fMassVsPtVsYME->Fill(mass,pt,rapid);
 903     }
 904   }
 905   return;
 906 }
 907 //________________________________________________________________________
 908 void AliAnalysisTaskCombinHF::FillMEHistosLS(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau, Int_t charge){
 909   // Fill histos for candidates in MixedEvents
 910
 911   tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
 912   Double_t pt = tmpRD->Pt();
 913   Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
 914   Double_t mass=TMath::Sqrt(minv2);
 915
 916   if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
 917     Double_t rapid = tmpRD->Y(pdgD);
 918     if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
 919       if(charge>0) fMassVsPtVsYMELSpp->Fill(mass,pt,rapid);
 920       if(charge<0) fMassVsPtVsYMELSmm->Fill(mass,pt,rapid);
 921     }
 922   }
 923   return;
 924 }
 925 //________________________________________________________________________
 926 Bool_t AliAnalysisTaskCombinHF::IsTrackSelected(AliAODTrack* track){
 927   // track selection cuts
 928
 929   if(track->Charge()==0) return kFALSE;
 930   if(track->GetID()<0&&!fKeepNegID)return kFALSE;
 931   if(!(track->TestFilterMask(fFilterMask))) return kFALSE;
 932   if(!SelectAODTrack(track,fTrackCutsAll)) return kFALSE;
 933   return kTRUE;
 934 }
 935
 936 //________________________________________________________________________
 937 Bool_t AliAnalysisTaskCombinHF::IsKaon(AliAODTrack* track){
 938   // kaon selection cuts
 939
 940   if(!fPidHF) return kTRUE;
 941   Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);
 942   Double_t mom=track->P();
 943   if(SelectAODTrack(track,fTrackCutsKaon)) {
 944     if(isKaon>=1)    return kTRUE;
 945     if(isKaon<=-1) return kFALSE;
 946     switch(fPIDselCaseZero){// isKaon=0
 947       case 0:
 948       {
 949         return kTRUE;// always accept
 950       }
 951         break;
 952       case 1:
 953       {
 954         if(isKaon>=0 && track->P()>fmaxPforIDKaon)return kTRUE;// accept only if in a compatibility band starting from p=fmaxPforIDKaon
 955       }
 956         break;
 957       case 2:
 958       {
 959         if(track->P()>fmaxPforIDKaon)return kTRUE;
 960         AliPIDResponse *pidResp=fPidHF->GetPidResponse();// more elaborated strategy: asymmetric cuts, with fix momenta and nsigma ranges for the moment
 961         Double_t nsigma=pidResp->NumberOfSigmasTPC(track,AliPID::kKaon);
 962         if(nsigma>-2.&& nsigma<3. && mom<0.6)isKaon=1;
 963         else if(nsigma>-1.&& nsigma<3.&& mom<0.8)isKaon=1;
 964         if(isKaon==1)return kTRUE;
 965       }
 966         break;
 967       default:
 968       {
 969         AliWarning(Form("WRONG CASE OF PID STRATEGY SELECTED: %d (can range from 0 to 2)",fPIDselCaseZero));
 970         return kFALSE;// actually case 0 could be set as the default and return kTRUE
 971       }
 972     }
 973   }
 974
 975   return kFALSE;
 976 }
 977 //_______________________________________________________________________
 978 Bool_t AliAnalysisTaskCombinHF::IsPion(AliAODTrack* track){
 979   // pion selection cuts
 980
 981   if(!fPidHF) return kTRUE;
 982   Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
 983   Double_t mom=track->P();
 984   if(SelectAODTrack(track,fTrackCutsPion)) {
 985     if(isPion>=1)    return kTRUE;
 986     if(isPion<=-1) return kFALSE;
 987     switch(fPIDselCaseZero){// isPion=0
 988       case 0:
 989       {
 990         return kTRUE;// always accept
 991       }
 992         break;
 993       case 1:
 994       {
 995         if(track->P()>fmaxPforIDPion)return kTRUE;// accept only if in a compatibility band starting from p=fmaxPforIDPion
 996       }
 997         break;
 998       case 2:
 999       {
1000         // more elaborated strategy: asymmetric cuts, with fix momenta and nsigma ranges for the moment
1001         if(track->P()>fmaxPforIDPion)return kTRUE;
1002         AliPIDResponse *pidResp=fPidHF->GetPidResponse();
1003         Double_t nsigma=pidResp->NumberOfSigmasTPC(track,AliPID::kPion);
1004         if(nsigma<2.&& nsigma>-3. && mom<0.6)isPion=1;
1005         else if(nsigma<1. && nsigma> -3. && mom<0.8)isPion=1;
1006         if(isPion==1)return kTRUE;
1007       }
1008         break;
1009       default:
1010       {
1011         AliWarning(Form("WRONG CASE OF PID STRATEGY SELECTED: %d (can range from 0 to 2)",fPIDselCaseZero));
1012         return kFALSE;// actually case 0 could be set as the default and return kTRUE
1013       }
1014     }
1015   }
1016
1017   return kFALSE;
1018 }
1019
1020 //________________________________________________________________________
1021 Bool_t AliAnalysisTaskCombinHF::SelectAODTrack(AliAODTrack *track, AliESDtrackCuts *cuts){
1022   // AOD track selection
1023
1024   if(!cuts) return kTRUE;
1025
1026   AliESDtrack esdTrack(track);
1027   // set the TPC cluster info
1028   esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
1029   esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
1030   esdTrack.SetTPCPointsF(track->GetTPCNclsF());
1031   if(!cuts->IsSelected(&esdTrack)) return kFALSE;
1032
1033   return kTRUE;
1034 }
1035
1036 //_________________________________________________________________
1037 Bool_t AliAnalysisTaskCombinHF::CheckAcceptance(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
1038   // check if the decay products are in the good eta and pt range
1039   for (Int_t iProng = 0; iProng<nProng; iProng++){
1040     AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
1041     if(!mcPartDaughter) return kFALSE;
1042     Double_t eta = mcPartDaughter->Eta();
1043     Double_t pt = mcPartDaughter->Pt();
1044     if (TMath::Abs(eta) > fEtaAccCut || pt < fPtAccCut) return kFALSE;
1045   }
1046   return kTRUE;
1047 }
1048 //_________________________________________________________________
1049 Int_t AliAnalysisTaskCombinHF::GetPoolIndex(Double_t zvert, Double_t mult){
1050   // check in which of the pools the current event falls
1051   if(!fzVertPoolLims || !fMultPoolLims) return 0;
1052   Int_t theBinZ=TMath::BinarySearch(fNzVertPoolsLimSize,fzVertPoolLims,zvert);
1053   if(theBinZ<0 || theBinZ>=fNzVertPoolsLimSize) return -1;
1054   Int_t theBinM=TMath::BinarySearch(fNMultPoolsLimSize,fMultPoolLims,mult);
1055   if(theBinM<0 || theBinM>=fNMultPoolsLimSize) return -1;
1056   return fNMultPools*theBinZ+theBinM;
1057 }
1058 //_________________________________________________________________
1059 void AliAnalysisTaskCombinHF::ResetPool(Int_t poolIndex){
1060   // delete the contets of the pool
1061   if(poolIndex<0 || poolIndex>=fNOfPools) return;
1062   delete fEventBuffer[poolIndex];
1063   fEventBuffer[poolIndex]=new TTree(Form("EventBuffer_%d",poolIndex), "Temporary buffer for event mixing");
1064   fEventBuffer[poolIndex]->Branch("zVertex", &fVtxZ);
1065   fEventBuffer[poolIndex]->Branch("multiplicity", &fMultiplicity);
1066   fEventBuffer[poolIndex]->Branch("eventInfo", "TObjString",&fEventInfo);
1067   fEventBuffer[poolIndex]->Branch("karray", "TObjArray", &fKaonTracks);
1068   fEventBuffer[poolIndex]->Branch("parray", "TObjArray", &fPionTracks);
1069   return;
1070 }
1071 //_________________________________________________________________
1072 Bool_t AliAnalysisTaskCombinHF::CanBeMixed(Double_t zv1, Double_t zv2, Double_t mult1, Double_t mult2){
1073   // check mixing
1074   if(TMath::Abs(zv2-zv1)>fMaxzVertDistForMix) return kFALSE;
1075   if(TMath::Abs(mult2-mult1)>fMaxMultDiffForMix) return kFALSE;
1076   return kTRUE;
1077 }
1078 //_________________________________________________________________
1079 void AliAnalysisTaskCombinHF::DoMixingWithCuts(){
1080   // perform mixed event analysis
1081
1082   if(fDoEventMixing==0) return;
1083   Int_t nEvents=fEventBuffer[0]->GetEntries();
1084   if(fDebug > 1) printf("AnalysisTaskCombinHF::DoMixingWithCuts Start Event Mixing of %d events\n",nEvents);
1085
1086   TObjArray* karray=0x0;
1087   TObjArray* parray=0x0;
1088   Double_t zVertex,mult;
1089   TObjString* eventInfo=0x0;
1090   fEventBuffer[0]->SetBranchAddress("karray", &karray);
1091   fEventBuffer[0]->SetBranchAddress("parray", &parray);
1092   fEventBuffer[0]->SetBranchAddress("eventInfo",&eventInfo);
1093   fEventBuffer[0]->SetBranchAddress("zVertex", &zVertex);
1094   fEventBuffer[0]->SetBranchAddress("multiplicity", &mult);
1095   Double_t d02[2]={0.,0.};
1096   Double_t d03[3]={0.,0.,0.};
1097   AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
1098   AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
1099   UInt_t pdg0[2]={321,211};
1100   //  UInt_t pdgp[3]={321,211,211};
1101   Double_t px[3],py[3],pz[3];
1102   Int_t evId1,esdId1,nk1,np1;
1103   Int_t evId2,esdId2,nk2,np2;
1104
1105   for(Int_t iEv1=0; iEv1<nEvents; iEv1++){
1106     fEventBuffer[0]->GetEvent(iEv1);
1107     TObjArray* karray1=(TObjArray*)karray->Clone();
1108     Double_t zVertex1=zVertex;
1109     Double_t mult1=mult;
1110     Int_t nKaons=karray1->GetEntries();
1111     Int_t nPionsForCheck=parray->GetEntries();
1112     sscanf((eventInfo->String()).Data(),"Ev%d_esd%d_Pi%d_K%d",&evId1,&esdId1,&np1,&nk1);
1113     if(nk1!=nKaons || np1!=nPionsForCheck){
1114       printf("AnalysisTaskCombinHF::DoMixingWithCuts ERROR: read event does not match to the stored one\n");
1115       delete karray1;
1116       continue;
1117     }
1118     for(Int_t iEv2=0; iEv2<fNumberOfEventsForMixing; iEv2++){
1119       Int_t iToMix=iEv1+iEv2+1;
1120       if(iEv1>=(nEvents-fNumberOfEventsForMixing)) iToMix=iEv1-iEv2-1;
1121       if(iToMix<0) continue;
1122       if(iToMix==iEv1) continue;
1123       if(iToMix<iEv1) continue;
1124       fEventBuffer[0]->GetEvent(iToMix);
1125       Double_t zVertex2=zVertex;
1126       Double_t mult2=mult;
1127       if(TMath::Abs(zVertex2-zVertex1)<0.0001 && TMath::Abs(mult2-mult1)<0.001){
1128         printf("AnalysisTaskCombinHF::DoMixingWithCuts ERROR: same event in mixing??? %d %d   %f %f  %f %f\n",iEv1,iEv2,zVertex1,zVertex2,mult1,mult2);
1129         continue;
1130       }
1131       TObjArray* parray2=(TObjArray*)parray->Clone();
1132       Int_t nPions=parray2->GetEntries();
1133       Int_t nKaonsForCheck=karray->GetEntries();
1134       sscanf((eventInfo->String()).Data(),"Ev%d_esd%d_Pi%d_K%d",&evId2,&esdId2,&np2,&nk2);
1135       if(nk2!=nKaonsForCheck || np2!=nPions){
1136         printf("AnalysisTaskCombinHF::DoMixingWithCuts ERROR: read event does not match to the stored one\n");
1137         delete parray2;
1138         continue;
1139       }
1140       if(evId2==evId1 && esdId2==esdId1){
1141         printf("AnalysisTaskCombinHF::DoMixingWithCuts ERROR: same event in mixing??? %d %d   nK=%d %d  nPi=%d %d\n",evId1,evId2,nKaons,nKaonsForCheck,nPionsForCheck,nPions);
1142         delete parray2;
1143         continue;
1144       }
1145       if(CanBeMixed(zVertex1,zVertex2,mult1,mult2)){
1146         for(Int_t iTr1=0; iTr1<nKaons; iTr1++){
1147           TLorentzVector* trK=(TLorentzVector*)karray1->At(iTr1);
1148           Double_t chargeK=trK->T();
1149           px[0] = trK->Px();
1150           py[0] = trK->Py();
1151           pz[0] = trK->Pz();
1152           for(Int_t iTr2=0; iTr2<nPions; iTr2++){
1153             TLorentzVector* trPi1=(TLorentzVector*)parray2->At(iTr2);
1154             Double_t chargePi1=trPi1->T();
1155             px[1] = trPi1->Px();
1156             py[1] = trPi1->Py();
1157             pz[1] = trPi1->Pz();
1158             if(fMeson==kDzero && chargePi1*chargeK<0){
1159               FillMEHistos(421,2,tmpRD2,px,py,pz,pdg0);
1160             }
1161             if(fMeson==kDzero && chargePi1*chargeK>0){
1162               FillMEHistosLS(421,2,tmpRD2,px,py,pz,pdg0,(Int_t)chargePi1);
1163             }
1164           }
1165         }
1166       }
1167       delete parray2;
1168     }
1169     delete karray1;
1170   }
1171   delete tmpRD2;
1172   delete tmpRD3;
1173 }
1174 //_________________________________________________________________
1175 void AliAnalysisTaskCombinHF::DoMixingWithPools(Int_t poolIndex){
1176   // perform mixed event analysis
1177
1178   if(fDoEventMixing==0) return;
1179   if(poolIndex<0 || poolIndex>fNzVertPools*fNMultPools) return;
1180
1181   Int_t nEvents=fEventBuffer[poolIndex]->GetEntries();
1182   if(fDebug > 1) printf("AliAnalysisTaskCombinHF::DoMixingWithPools Start Event Mixing of %d events\n",nEvents);
1183   TObjArray* karray=0x0;
1184   TObjArray* parray=0x0;
1185   Double_t zVertex,mult;
1186   TObjString* eventInfo=0x0;
1187   fEventBuffer[poolIndex]->SetBranchAddress("karray", &karray);
1188   fEventBuffer[poolIndex]->SetBranchAddress("parray", &parray);
1189   fEventBuffer[poolIndex]->SetBranchAddress("eventInfo",&eventInfo);
1190   fEventBuffer[poolIndex]->SetBranchAddress("zVertex", &zVertex);
1191   fEventBuffer[poolIndex]->SetBranchAddress("multiplicity", &mult);
1192
1193   // dummy values of track impact parameter, needed to build an AliAODRecoDecay object
1194   Double_t d02[2]={0.,0.};
1195   Double_t d03[3]={0.,0.,0.};
1196   AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
1197   AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
1198   UInt_t pdg0[2]={321,211};
1199   UInt_t pdgp[3]={321,211,211};
1200   Double_t px[3],py[3],pz[3];
1201   Int_t evId1,esdId1,nk1,np1;
1202   Int_t evId2,esdId2,nk2,np2;
1203
1204   for(Int_t iEv1=0; iEv1<nEvents; iEv1++){
1205     fEventBuffer[poolIndex]->GetEvent(iEv1);
1206     TObjArray* karray1=(TObjArray*)karray->Clone();
1207     Double_t zVertex1=zVertex;
1208     Double_t mult1=mult;
1209     Int_t nKaons=karray1->GetEntries();
1210     Int_t nPionsForCheck=parray->GetEntries();
1211     sscanf((eventInfo->String()).Data(),"Ev%d_esd%d_Pi%d_K%d",&evId1,&esdId1,&np1,&nk1);
1212     if(nk1!=nKaons || np1!=nPionsForCheck){
1213       printf("AliAnalysisTaskCombinHF::DoMixingWithPools ERROR: read event does not match to the stored one\n");
1214       delete karray1;
1215       continue;
1216     }
1217     for(Int_t iEv2=0; iEv2<nEvents; iEv2++){
1218       if(iEv2==iEv1) continue;
1219       fEventBuffer[poolIndex]->GetEvent(iEv2);
1220       Double_t zVertex2=zVertex;
1221       Double_t mult2=mult;
1222       if(TMath::Abs(zVertex2-zVertex1)<0.0001 && TMath::Abs(mult2-mult1)<0.001){
1223         printf("AliAnalysisTaskCombinHF::DoMixingWithPools ERROR: same event in mixing??? %d %d   %f %f  %f %f\n",iEv1,iEv2,zVertex1,zVertex2,mult1,mult2);
1224         continue;
1225       }
1226       TObjArray* parray2=(TObjArray*)parray->Clone();
1227       Int_t nPions=parray2->GetEntries();
1228       Int_t nKaonsForCheck=karray->GetEntries();
1229       sscanf((eventInfo->String()).Data(),"Ev%d_esd%d_Pi%d_K%d",&evId2,&esdId2,&np2,&nk2);
1230       if(nk2!=nKaonsForCheck || np2!=nPions){
1231         printf("AliAnalysisTaskCombinHF::DoMixingWithPools ERROR: read event does not match to the stored one\n");
1232         delete parray2;
1233         continue;
1234       }
1235       if(evId2==evId1 && esdId2==esdId1){
1236         printf("AliAnalysisTaskCombinHF::DoMixingWithPools ERROR: same event in mixing??? %d %d   nK=%d %d  nPi=%d %d\n",evId1,evId2,nKaons,nKaonsForCheck,nPionsForCheck,nPions);
1237         delete parray2;
1238         continue;
1239       }
1240       TObjArray* parray3=0x0;
1241       Int_t nPions3=0;
1242       if(fMeson!=kDzero){
1243         Int_t iEv3=iEv2+1;
1244         if(iEv3==iEv1) iEv3=iEv2+2;
1245         if(iEv3>=nEvents) iEv3=iEv2-3;
1246         if(nEvents==2) iEv3=iEv1;
1247         if(iEv3<0) iEv3=iEv2-1;
1248         fEventBuffer[poolIndex]->GetEvent(iEv3);
1249         parray3=(TObjArray*)parray->Clone();
1250         nPions3=parray3->GetEntries();
1251       }
1252       for(Int_t iTr1=0; iTr1<nKaons; iTr1++){
1253         TLorentzVector* trK=(TLorentzVector*)karray1->At(iTr1);
1254         Double_t chargeK=trK->T();
1255         px[0] = trK->Px();
1256         py[0] = trK->Py();
1257         pz[0] = trK->Pz();
1258         for(Int_t iTr2=0; iTr2<nPions; iTr2++){
1259           TLorentzVector* trPi1=(TLorentzVector*)parray2->At(iTr2);
1260           Double_t chargePi1=trPi1->T();
1261           px[1] = trPi1->Px();
1262           py[1] = trPi1->Py();
1263           pz[1] = trPi1->Pz();
1264           if(chargePi1*chargeK<0){
1265             if(fMeson==kDzero){
1266               FillMEHistos(421,2,tmpRD2,px,py,pz,pdg0);
1267             }else{
1268               if(parray3){
1269                 for(Int_t iTr3=iTr2+1; iTr3<nPions3; iTr3++){
1270                   TLorentzVector* trPi2=(TLorentzVector*)parray3->At(iTr3);
1271                   Double_t chargePi2=trPi2->T();
1272                   px[2] = trPi2->Px();
1273                   py[2] = trPi2->Py();
1274                   pz[2] = trPi2->Pz();
1275                   if(chargePi2*chargeK<0){
1276                     FillMEHistos(411,3,tmpRD3,px,py,pz,pdgp);
1277                   }
1278                 }
1279               }
1280             }
1281           }else if(chargePi1*chargeK>0){
1282             if(fMeson==kDzero) FillMEHistosLS(421,2,tmpRD2,px,py,pz,pdg0,(Int_t)chargePi1);
1283           }
1284         }
1285       }
1286       delete parray3;
1287       delete parray2;
1288     }
1289     delete karray1;
1290   }
1291   delete tmpRD2;
1292   delete tmpRD3;
1293 }
1294 //_________________________________________________________________
1295 void AliAnalysisTaskCombinHF::FinishTaskOutput()
1296 {
1297   // perform mixed event analysis
1298   if(fDoEventMixing==0) return;
1299   printf("AliAnalysisTaskCombinHF: FinishTaskOutput\n");
1300
1301   if(fDoEventMixing==1){
1302     for(Int_t i=0; i<fNOfPools; i++){
1303       Int_t nEvents=fEventBuffer[i]->GetEntries();
1304       if(nEvents>1) DoMixingWithPools(i);
1305     }
1306   }else if(fDoEventMixing==2){
1307     DoMixingWithCuts();
1308   }
1309 }
1310 //_________________________________________________________________
1311 void AliAnalysisTaskCombinHF::Terminate(Option_t */*option*/)
1312 {
1313   // Terminate analysis
1314   //
1315   if(fDebug > 1) printf("AliAnalysisTaskCombinHF: Terminate() \n");
1316   fOutput = dynamic_cast<TList*> (GetOutputData(1));
1317   if (!fOutput) {
1318     printf("ERROR: fOutput not available\n");
1319     return;
1320   }
1321   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
1322   if(fHistNEvents){
1323     printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
1324   }else{
1325     printf("ERROR: fHistNEvents not available\n");
1326     return;
1327   }
1328   return;
1329 }
1330