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   fPtVsYGen(0x0),
  57   fPtVsYGenLargeAcc(0x0),
  58   fPtVsYGenLimAcc(0x0),
  59   fPtVsYGenAcc(0x0),
  60   fPtVsYReco(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   fFilterMask(BIT(4)),
  74   fTrackCutsAll(0x0),
  75   fTrackCutsPion(0x0),
  76   fTrackCutsKaon(0x0),
  77   fPidHF(new AliAODPidHF()),
  78   fAnalysisCuts(0x0),
  79   fMinMass(1.720),
  80   fMaxMass(2.150),
  81   fMaxPt(10.),
  82   fPtBinWidth(0.5),
  83   fEtaAccCut(0.9),
  84   fPtAccCut(0.1),
  85   fNRotations(9),
  86   fMinAngleForRot(5*TMath::Pi()/6),
  87   fMaxAngleForRot(7*TMath::Pi()/6),
  88   fMinAngleForRot3(2*TMath::Pi()/6),
  89   fMaxAngleForRot3(4*TMath::Pi()/6),
  90   fCounter(0x0),
  91   fMeson(kDzero),
  92   fReadMC(kFALSE),
  93   fPromptFeeddown(kPrompt),
  94   fGoUpToQuark(kTRUE),
  95   fFullAnalysis(0),
  96   fPIDstrategy(knSigma),
  97   fmaxPforIDPion(0.8),
  98   fmaxPforIDKaon(2.),
  99   fKeepNegID(kFALSE),
 100   fPIDselCaseZero(0),
 101   fBayesThresKaon(0.4),
 102   fBayesThresPion(0.4),
 103   fDoEventMixing(kTRUE),
 104   fMaxNumberOfEventsForMixing(20),
 105   fMaxzVertDistForMix(20.),
 106   fMaxMultDiffForMix(9999.),
 107   fUsePoolsZ(kFALSE),
 108   fUsePoolsM(kFALSE),
 109   fNzVertPools(1),
 110   fNzVertPoolsLimSize(2),
 111   fzVertPoolLims(0x0),
 112   fNMultPools(1),
 113   fNMultPoolsLimSize(2),
 114   fMultPoolLims(0x0),
 115   fEventBuffer(0x0),
 116   fVtxZ(0),
 117   fMultiplicity(0),
 118   fKaonTracks(0x0),
 119   fPionTracks(0x0)
 120 {
 121   // default constructor
 122 }
 123
 124 //________________________________________________________________________
 125 AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF(Int_t meson, AliRDHFCuts* analysiscuts):
 126   AliAnalysisTaskSE("DmesonCombin"),
 127   fOutput(0x0),
 128   fHistNEvents(0x0),
 129   fHistTrackStatus(0x0),
 130   fHistCheckOrigin(0x0),
 131   fHistCheckOriginSel(0x0),
 132   fHistCheckDecChan(0x0),
 133   fHistCheckDecChanAcc(0x0),
 134   fPtVsYGen(0x0),
 135   fPtVsYGenLargeAcc(0x0),
 136   fPtVsYGenLimAcc(0x0),
 137   fPtVsYGenAcc(0x0),
 138   fPtVsYReco(0x0),
 139   fMassVsPtVsY(0x0),
 140   fMassVsPtVsYRot(0x0),
 141   fMassVsPtVsYLSpp(0x0),
 142   fMassVsPtVsYLSmm(0x0),
 143   fMassVsPtVsYSig(0x0),
 144   fMassVsPtVsYRefl(0x0),
 145   fMassVsPtVsYBkg(0x0),
 146   fNSelected(0x0),
 147   fNormRotated(0x0),
 148   fDeltaMass(0x0),
 149   fDeltaMassFullAnalysis(0x0),
 150   fMassVsPtVsYME(0x0),
 151   fFilterMask(BIT(4)),
 152   fTrackCutsAll(0x0),
 153   fTrackCutsPion(0x0),
 154   fTrackCutsKaon(0x0),
 155   fPidHF(new AliAODPidHF()),
 156   fAnalysisCuts(analysiscuts),
 157   fMinMass(1.720),
 158   fMaxMass(2.150),
 159   fMaxPt(10.),
 160   fPtBinWidth(0.5),
 161   fEtaAccCut(0.9),
 162   fPtAccCut(0.1),
 163   fNRotations(9),
 164   fMinAngleForRot(5*TMath::Pi()/6),
 165   fMaxAngleForRot(7*TMath::Pi()/6),
 166   fMinAngleForRot3(2*TMath::Pi()/6),
 167   fMaxAngleForRot3(4*TMath::Pi()/6),
 168   fCounter(0x0),
 169   fMeson(meson),
 170   fReadMC(kFALSE),
 171   fPromptFeeddown(1),
 172   fGoUpToQuark(kTRUE),
 173   fFullAnalysis(0),
 174   fPIDstrategy(knSigma),
 175   fmaxPforIDPion(0.8),
 176   fmaxPforIDKaon(2.),
 177   fKeepNegID(kFALSE),
 178   fPIDselCaseZero(0),
 179   fBayesThresKaon(0.4),
 180   fBayesThresPion(0.4),
 181   fDoEventMixing(kTRUE),
 182   fMaxNumberOfEventsForMixing(20),
 183   fMaxzVertDistForMix(20.),
 184   fMaxMultDiffForMix(9999.),
 185   fUsePoolsZ(kFALSE),
 186   fUsePoolsM(kFALSE),
 187   fNzVertPools(1),
 188   fNzVertPoolsLimSize(2),
 189   fzVertPoolLims(0x0),
 190   fNMultPools(1),
 191   fNMultPoolsLimSize(2),
 192   fMultPoolLims(0x0),
 193   fEventBuffer(0x0),
 194   fVtxZ(0),
 195   fMultiplicity(0),
 196   fKaonTracks(0x0),
 197   fPionTracks(0x0)
 198 {
 199   // standard constructor
 200
 201   DefineOutput(1,TList::Class());  //My private output
 202   DefineOutput(2,AliNormalizationCounter::Class());
 203 }
 204
 205 //________________________________________________________________________
 206 AliAnalysisTaskCombinHF::~AliAnalysisTaskCombinHF()
 207 {
 208   //
 209   // Destructor
 210   //
 211   if(fOutput && !fOutput->IsOwner()){
 212     delete fHistNEvents;
 213     delete fHistTrackStatus;
 214     delete fHistCheckOrigin;
 215     delete fHistCheckOriginSel;
 216     delete fHistCheckDecChan;
 217     delete fHistCheckDecChanAcc;
 218     delete fPtVsYGen;
 219     delete fPtVsYGenLargeAcc;
 220     delete fPtVsYGenLimAcc;
 221     delete fPtVsYGenAcc;
 222     delete fPtVsYReco;
 223     delete fMassVsPtVsY;
 224     delete fMassVsPtVsYLSpp;
 225     delete fMassVsPtVsYLSmm;
 226     delete fMassVsPtVsYRot;
 227     delete fMassVsPtVsYSig;
 228     delete fMassVsPtVsYRefl;
 229     delete fMassVsPtVsYBkg;
 230     delete fNSelected;
 231     delete fNormRotated;
 232     delete fDeltaMass;
 233     delete fDeltaMassFullAnalysis;
 234     delete fMassVsPtVsYME;
 235   }
 236
 237   delete fOutput;
 238   delete fCounter;
 239   delete fTrackCutsAll;
 240   delete fTrackCutsPion;
 241   delete fTrackCutsKaon;
 242   delete fPidHF;
 243   delete fAnalysisCuts;
 244   if(fKaonTracks) fKaonTracks->Delete();
 245   if(fPionTracks) fPionTracks->Delete();
 246   delete fKaonTracks;
 247   delete fPionTracks;
 248   delete fEventBuffer;
 249   delete [] fzVertPoolLims;
 250   delete [] fMultPoolLims;
 251 }
 252
 253 //________________________________________________________________________
 254 void AliAnalysisTaskCombinHF::ConfigureZVertPools(Int_t nPools, Double_t*  zVertLimits)
 255 {
 256   // sets the pools for event mizing in zvertex
 257   if(fzVertPoolLims) delete [] fzVertPoolLims;
 258   fNzVertPools=nPools;
 259   fNzVertPoolsLimSize=nPools+1;
 260   fzVertPoolLims = new Double_t[fNzVertPoolsLimSize];
 261   for(Int_t ib=0; ib<fNzVertPoolsLimSize; ib++) fzVertPoolLims[ib]=zVertLimits[ib];
 262   fUsePoolsZ=kTRUE;
 263   return;
 264 }
 265 //________________________________________________________________________
 266 void AliAnalysisTaskCombinHF::ConfigureMultiplicityPools(Int_t nPools, Double_t*  multLimits)
 267 {
 268   // sets the pools for event mizing in zvertex
 269   if(fMultPoolLims) delete [] fMultPoolLims;
 270   fNMultPools=nPools;
 271   fNMultPoolsLimSize=nPools+1;
 272   fMultPoolLims = new Double_t[fNMultPoolsLimSize];
 273   for(Int_t ib=0; ib<nPools+1; ib++) fMultPoolLims[ib]=multLimits[ib];
 274   fUsePoolsM=kTRUE;
 275   return;
 276 }
 277 //________________________________________________________________________
 278 void AliAnalysisTaskCombinHF::UserCreateOutputObjects()
 279 {
 280   // Create the output container
 281   //
 282   if(fDebug > 1) printf("AnalysisTaskCombinHF::UserCreateOutputObjects() \n");
 283
 284   fOutput = new TList();
 285   fOutput->SetOwner();
 286   fOutput->SetName("OutputHistos");
 287
 288   fHistNEvents = new TH1F("hNEvents", "number of events ",8,-0.5,7.5);
 289   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 290   fHistNEvents->GetXaxis()->SetBinLabel(2,"n. passing IsEvSelected");
 291   fHistNEvents->GetXaxis()->SetBinLabel(3,"n. rejected due to trigger");
 292   fHistNEvents->GetXaxis()->SetBinLabel(4,"n. rejected due to not reco vertex");
 293   fHistNEvents->GetXaxis()->SetBinLabel(5,"n. rejected for contr vertex");
 294   fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected for vertex out of accept");
 295   fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected for pileup events");
 296   fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events");
 297
 298   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
 299   fHistNEvents->Sumw2();
 300   fHistNEvents->SetMinimum(0);
 301   fOutput->Add(fHistNEvents);
 302
 303   fHistTrackStatus  = new TH1F("hTrackStatus", "",8,-0.5,7.5);
 304   fHistTrackStatus->GetXaxis()->SetBinLabel(1,"Not OK");
 305   fHistTrackStatus->GetXaxis()->SetBinLabel(2,"Track OK");
 306   fHistTrackStatus->GetXaxis()->SetBinLabel(3,"Kaon, Not OK");
 307   fHistTrackStatus->GetXaxis()->SetBinLabel(4,"Kaon OK");
 308   fHistTrackStatus->GetXaxis()->SetBinLabel(5,"Pion, Not OK");
 309   fHistTrackStatus->GetXaxis()->SetBinLabel(6,"Pion OK");
 310   fHistTrackStatus->GetXaxis()->SetBinLabel(7,"Kaon||Pion, Not OK");
 311   fHistTrackStatus->GetXaxis()->SetBinLabel(8,"Kaon||Pion OK");
 312   fHistTrackStatus->GetXaxis()->SetNdivisions(1,kFALSE);
 313   fHistTrackStatus->Sumw2();
 314   fHistTrackStatus->SetMinimum(0);
 315   fOutput->Add(fHistTrackStatus);
 316
 317   Int_t nPtBins = (Int_t)(fMaxPt/fPtBinWidth+0.001);
 318   Double_t maxPt=fPtBinWidth*nPtBins;
 319
 320   if(fReadMC){
 321
 322     fHistCheckOrigin=new TH1F("hCheckOrigin","",7,-1.5,5.5);
 323     fHistCheckOrigin->Sumw2();
 324     fHistCheckOrigin->SetMinimum(0);
 325     fOutput->Add(fHistCheckOrigin);
 326
 327     fHistCheckOriginSel=new TH1F("hCheckOriginSel","",7,-1.5,5.5);
 328     fHistCheckOriginSel->Sumw2();
 329     fHistCheckOriginSel->SetMinimum(0);
 330     fOutput->Add(fHistCheckOriginSel);
 331
 332     fHistCheckDecChan=new TH1F("hCheckDecChan","",7,-2.5,4.5);
 333     fHistCheckDecChan->Sumw2();
 334     fHistCheckDecChan->SetMinimum(0);
 335     fOutput->Add(fHistCheckDecChan);
 336
 337     fHistCheckDecChanAcc=new TH1F("hCheckDecChanAcc","",7,-2.5,4.5);
 338     fHistCheckDecChanAcc->Sumw2();
 339     fHistCheckDecChanAcc->SetMinimum(0);
 340     fOutput->Add(fHistCheckDecChanAcc);
 341
 342     fPtVsYGen= new TH2F("hPtVsYGen","",nPtBins,0.,maxPt,20,-1.,1.);
 343     fPtVsYGen->Sumw2();
 344     fPtVsYGen->SetMinimum(0);
 345     fOutput->Add(fPtVsYGen);
 346
 347     fPtVsYGenLargeAcc= new TH2F("hPtVsYGenLargeAcc","",nPtBins,0.,maxPt,20,-1.,1.);
 348     fPtVsYGenLargeAcc->Sumw2();
 349     fPtVsYGenLargeAcc->SetMinimum(0);
 350     fOutput->Add(fPtVsYGenLargeAcc);
 351
 352     fPtVsYGenLimAcc= new TH2F("hPtVsYGenLimAcc","",nPtBins,0.,maxPt,20,-1.,1.);
 353     fPtVsYGenLimAcc->Sumw2();
 354     fPtVsYGenLimAcc->SetMinimum(0);
 355     fOutput->Add(fPtVsYGenLimAcc);
 356
 357     fPtVsYGenAcc= new TH2F("hPtVsYGenAcc","",nPtBins,0.,maxPt,20,-1.,1.);
 358     fPtVsYGenAcc->Sumw2();
 359     fPtVsYGenAcc->SetMinimum(0);
 360     fOutput->Add(fPtVsYGenAcc);
 361
 362     fPtVsYReco= new TH2F("hPtVsYReco","",nPtBins,0.,maxPt,20,-1.,1.);
 363     fPtVsYReco->Sumw2();
 364     fPtVsYReco->SetMinimum(0);
 365     fOutput->Add(fPtVsYReco);
 366   }
 367
 368
 369   Int_t nMassBins=fMaxMass*1000.-fMinMass*1000.;
 370   Double_t maxm=fMinMass+nMassBins*0.001;
 371   fMassVsPtVsY=new TH3F("hMassVsPtVsY","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
 372   fMassVsPtVsY->Sumw2();
 373   fMassVsPtVsY->SetMinimum(0);
 374   fOutput->Add(fMassVsPtVsY);
 375
 376   fMassVsPtVsYRot=new TH3F("hMassVsPtVsYRot","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
 377   fMassVsPtVsYRot->Sumw2();
 378   fMassVsPtVsYRot->SetMinimum(0);
 379   fOutput->Add(fMassVsPtVsYRot);
 380
 381   if(fMeson==kDzero){
 382     fMassVsPtVsYLSpp=new TH3F("hMassVsPtVsYLSpp","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
 383     fMassVsPtVsYLSpp->Sumw2();
 384     fMassVsPtVsYLSpp->SetMinimum(0);
 385     fOutput->Add(fMassVsPtVsYLSpp);
 386     fMassVsPtVsYLSmm=new TH3F("hMassVsPtVsYLSmm","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
 387     fMassVsPtVsYLSmm->Sumw2();
 388     fMassVsPtVsYLSmm->SetMinimum(0);
 389     fOutput->Add(fMassVsPtVsYLSmm);
 390   }
 391
 392   fMassVsPtVsYSig=new TH3F("hMassVsPtVsYSig","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
 393   fMassVsPtVsYSig->Sumw2();
 394   fMassVsPtVsYSig->SetMinimum(0);
 395   fOutput->Add(fMassVsPtVsYSig);
 396
 397   fMassVsPtVsYRefl=new TH3F("hMassVsPtVsYRefl","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
 398   fMassVsPtVsYRefl->Sumw2();
 399   fMassVsPtVsYRefl->SetMinimum(0);
 400   fOutput->Add(fMassVsPtVsYRefl);
 401
 402   fMassVsPtVsYBkg=new TH3F("hMassVsPtVsYBkg","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
 403   fMassVsPtVsYBkg->Sumw2();
 404   fMassVsPtVsYBkg->SetMinimum(0);
 405   fOutput->Add(fMassVsPtVsYBkg);
 406
 407   fNSelected=new TH1F("hNSelected","",100,-0.5,99.5);
 408   fNSelected->Sumw2();
 409   fNSelected->SetMinimum(0);
 410   fOutput->Add(fNSelected);
 411
 412   fNormRotated=new TH1F("hNormRotated","",11,-0.5,10.5);
 413   fNormRotated->Sumw2();
 414   fNormRotated->SetMinimum(0);
 415   fOutput->Add(fNormRotated);
 416
 417   fDeltaMass=new TH1F("hDeltaMass","",100,-0.4,0.4);
 418   fDeltaMass->Sumw2();
 419   fDeltaMass->SetMinimum(0);
 420   fOutput->Add(fDeltaMass);
 421
 422   Int_t binSparseDMassRot[5]={nMassBins,100,24,40,20};
 423   Double_t edgeLowSparseDMassRot[5]={fMinMass,-0.4,0.,-4.,0};
 424   Double_t edgeHighSparseDMassRot[5]={maxm,0.4,12.,4.,3.14};
 425   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);
 426   fOutput->Add(fDeltaMassFullAnalysis);
 427
 428   fMassVsPtVsYME=new TH3F("hMassVsPtVsYME","",nMassBins,fMinMass,maxm,nPtBins,0.,maxPt,20,-1.,1.);
 429   fMassVsPtVsYME->Sumw2();
 430   fMassVsPtVsYME->SetMinimum(0);
 431   fOutput->Add(fMassVsPtVsYME);
 432
 433
 434   //Counter for Normalization
 435   fCounter = new AliNormalizationCounter("NormalizationCounter");
 436   fCounter->Init();
 437
 438   fKaonTracks = new TObjArray();
 439   fPionTracks=new TObjArray();
 440   fKaonTracks->SetOwner();
 441   fPionTracks->SetOwner();
 442   fEventBuffer = new TTree("EventBuffer", "Temporary buffer for event mixing");
 443   fEventBuffer->Branch("zVertex", &fVtxZ);
 444   fEventBuffer->Branch("multiplicity", &fMultiplicity);
 445   fEventBuffer->Branch("karray", "TObjArray", &fKaonTracks);
 446   fEventBuffer->Branch("parray", "TObjArray", &fPionTracks);
 447
 448   PostData(1,fOutput);
 449   PostData(2,fCounter);
 450 }
 451
 452 //________________________________________________________________________
 453 void AliAnalysisTaskCombinHF::UserExec(Option_t */*option*/){
 454   //Build the 3-track combinatorics (+-+ and -+-) for D+->Kpipi decays
 455
 456   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 457   if(!aod && AODEvent() && IsStandardAOD()) {
 458     // In case there is an AOD handler writing a standard AOD, use the AOD
 459     // event in memory rather than the input (ESD) event.
 460     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 461   }
 462   if(!aod){
 463     printf("AliAnalysisTaskCombinHF::UserExec: AOD not found!\n");
 464     return;
 465   }
 466
 467   // fix for temporary bug in ESDfilter
 468   // the AODs with null vertex pointer didn't pass the PhysSel
 469   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
 470   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 471   AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
 472   AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
 473   fPidHF->SetPidResponse(pidResp);
 474
 475
 476   fHistNEvents->Fill(0); // count event
 477   // Post the data already here
 478   PostData(1,fOutput);
 479
 480   fCounter->StoreEvent(aod,fAnalysisCuts,fReadMC);
 481
 482   Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
 483   if(fAnalysisCuts->IsEventRejectedDueToTrigger())fHistNEvents->Fill(2);
 484   if(fAnalysisCuts->IsEventRejectedDueToNotRecoVertex())fHistNEvents->Fill(3);
 485   if(fAnalysisCuts->IsEventRejectedDueToVertexContributors())fHistNEvents->Fill(4);
 486   if(fAnalysisCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion())fHistNEvents->Fill(5);
 487   //  if(fAnalysisCuts->IsEventRejectedDueToPileup())fHistNEvents->Fill(6);
 488   if(fAnalysisCuts->IsEventRejectedDueToCentrality()) fHistNEvents->Fill(7);
 489
 490   if(!isEvSel)return;
 491
 492   fHistNEvents->Fill(1);
 493
 494   TClonesArray *arrayMC=0;
 495   AliAODMCHeader *mcHeader=0;
 496   if(fReadMC){
 497     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 498     if(!arrayMC) {
 499       printf("AliAnalysisTaskCombinHF::UserExec: MC particles branch not found!\n");
 500       return;
 501     }
 502
 503     // load MC header
 504     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 505     if(!mcHeader) {
 506       printf("AliAnalysisTaskCombinHF::UserExec: MC header branch not found!\n");
 507       return;
 508     }
 509     FillGenHistos(arrayMC);
 510   }
 511
 512
 513   Int_t ntracks=aod->GetNTracks();
 514   fVtxZ = aod->GetPrimaryVertex()->GetZ();
 515   fMultiplicity = AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
 516
 517   // select and flag tracks
 518   UChar_t* status = new UChar_t[ntracks];
 519   for(Int_t iTr=0; iTr<ntracks; iTr++){
 520     status[iTr]=0;
 521     AliAODTrack* track=aod->GetTrack(iTr);
 522     if(IsTrackSelected(track)) status[iTr]+=1;
 523
 524     // PID
 525     if (fPIDstrategy == knSigma) {
 526       // nsigma PID
 527       if(IsKaon(track)) status[iTr]+=2;
 528       if(IsPion(track)) status[iTr]+=4;
 529     }
 530     else if (fPIDstrategy == kBayesianMaxProb || fPIDstrategy == kBayesianThres) {
 531       // Bayesian PID
 532       Double_t *weights = new Double_t[AliPID::kSPECIES];
 533       fPidHF->GetPidCombined()->ComputeProbabilities(track, fPidHF->GetPidResponse(), weights);
 534       if (fPIDstrategy == kBayesianMaxProb) {
 535         if (TMath::MaxElement(AliPID::kSPECIES, weights) == weights[AliPID::kKaon]) status[iTr] += 2;
 536         if (TMath::MaxElement(AliPID::kSPECIES, weights) == weights[AliPID::kPion]) status[iTr] += 4;
 537       }
 538       if (fPIDstrategy == kBayesianThres) {
 539         if (weights[AliPID::kKaon] > fBayesThresKaon) status[iTr] += 2;
 540         if (weights[AliPID::kPion] > fBayesThresPion) status[iTr] += 4;
 541       }
 542       delete[] weights;
 543     }
 544
 545     fHistTrackStatus->Fill(status[iTr]);
 546   }
 547
 548   // build the combinatorics
 549   Int_t nSelected=0;
 550   Int_t nFiltered=0;
 551   Double_t dummypos[3]={0.,0.,0.};
 552   AliAODVertex* v2=new AliAODVertex(dummypos,999.,-1,2);
 553   AliAODVertex* v3=new AliAODVertex(dummypos,999.,-1,3);
 554   // dummy values of track impact parameter, needed to build an AliAODRecoDecay object
 555   Double_t d02[2]={0.,0.};
 556   Double_t d03[3]={0.,0.,0.};
 557   AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
 558   AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
 559   UInt_t pdg0[2]={321,211};
 560   UInt_t pdgp[3]={321,211,211};
 561   //  UInt_t pdgs[3]={321,321,211};
 562   Double_t tmpp[3];
 563   Double_t px[3],py[3],pz[3];
 564   Int_t dgLabels[3];
 565   fKaonTracks->Delete();
 566   fPionTracks->Delete();
 567
 568   for(Int_t iTr1=0; iTr1<ntracks; iTr1++){
 569     AliAODTrack* trK=aod->GetTrack(iTr1);
 570     if((status[iTr1] & 1)==0) continue;
 571     if(fDoEventMixing){
 572       if(status[iTr1] & 2) fKaonTracks->AddLast(new TLorentzVector(trK->Px(),trK->Py(),trK->Pz(),trK->Charge()));
 573       if(status[iTr1] & 4) fPionTracks->AddLast(new TLorentzVector(trK->Px(),trK->Py(),trK->Pz(),trK->Charge()));
 574     }
 575     if((status[iTr1] & 2)==0) continue;
 576     Int_t chargeK=trK->Charge();
 577     trK->GetPxPyPz(tmpp);
 578     px[0] = tmpp[0];
 579     py[0] = tmpp[1];
 580     pz[0] = tmpp[2];
 581     dgLabels[0]=trK->GetLabel();
 582     for(Int_t iTr2=0; iTr2<ntracks; iTr2++){
 583       if((status[iTr2] & 1)==0) continue;
 584       if((status[iTr2] & 4)==0) continue;
 585       if(iTr1==iTr2) continue;
 586       AliAODTrack* trPi1=aod->GetTrack(iTr2);
 587       Int_t chargePi1=trPi1->Charge();
 588       trPi1->GetPxPyPz(tmpp);
 589       px[1] = tmpp[0];
 590       py[1] = tmpp[1];
 591       pz[1] = tmpp[2];
 592       dgLabels[1]=trPi1->GetLabel();
 593       if(chargePi1==chargeK){
 594         if(fMeson==kDzero) FillLSHistos(421,2,tmpRD2,px,py,pz,pdg0,chargePi1);
 595         continue;
 596       }
 597       if(fMeson==kDzero){
 598         nFiltered++;
 599         v2->AddDaughter(trK);
 600         v2->AddDaughter(trPi1);
 601         tmpRD2->SetSecondaryVtx(v2);
 602         Bool_t ok=FillHistos(421,2,tmpRD2,px,py,pz,pdg0,arrayMC,dgLabels);
 603         v2->RemoveDaughters();
 604         if(ok) nSelected++;
 605       }else{
 606         for(Int_t iTr3=iTr2+1; iTr3<ntracks; iTr3++){
 607           if((status[iTr3] & 1)==0) continue;
 608           if((status[iTr3] & 4)==0) continue;
 609           if(iTr1==iTr3) continue;
 610           AliAODTrack* trPi2=aod->GetTrack(iTr3);
 611           Int_t chargePi2=trPi2->Charge();
 612           if(chargePi2==chargeK) continue;
 613           nFiltered++;
 614           trPi2->GetPxPyPz(tmpp);
 615           px[2] = tmpp[0];
 616           py[2] = tmpp[1];
 617           pz[2] = tmpp[2];
 618           dgLabels[2]=trPi2->GetLabel();
 619           v3->AddDaughter(trK);
 620           v3->AddDaughter(trPi1);
 621           v3->AddDaughter(trPi2);
 622           tmpRD3->SetSecondaryVtx(v3);
 623           Bool_t ok=FillHistos(411,3,tmpRD3,px,py,pz,pdgp,arrayMC,dgLabels);
 624           v3->RemoveDaughters();
 625           if(ok) nSelected++;
 626         }
 627       }
 628     }
 629   }
 630
 631   delete [] status;
 632   delete v2;
 633   delete v3;
 634   delete tmpRD2;
 635   delete tmpRD3;
 636
 637   fNSelected->Fill(nSelected);
 638
 639   fCounter->StoreCandidates(aod,nFiltered,kTRUE);
 640   fCounter->StoreCandidates(aod,nSelected,kFALSE);
 641
 642   if(fDoEventMixing) fEventBuffer->Fill();
 643   PostData(1,fOutput);
 644   PostData(2,fCounter);
 645
 646   return;
 647 }
 648
 649 //________________________________________________________________________
 650 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){
 651   // Fill histos for LS candidates
 652
 653   tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
 654   Double_t pt = tmpRD->Pt();
 655   Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
 656   if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
 657     Double_t rapid = tmpRD->Y(pdgD);
 658     if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
 659       if(charge>0) fMassVsPtVsYLSpp->Fill(TMath::Sqrt(minv2),pt,rapid);
 660       else fMassVsPtVsYLSmm->Fill(TMath::Sqrt(minv2),pt,rapid);
 661     }
 662   }
 663   return;
 664 }
 665
 666 //________________________________________________________________________
 667 void AliAnalysisTaskCombinHF::FillGenHistos(TClonesArray* arrayMC){
 668   // Fill histos with generated quantities
 669   Int_t totPart=arrayMC->GetEntriesFast();
 670   Int_t thePDG=411;
 671   Int_t nProng=3;
 672   if(fMeson==kDzero){
 673     thePDG=421;
 674     nProng=2;
 675   }
 676   for(Int_t ip=0; ip<totPart; ip++){
 677     AliAODMCParticle *part = (AliAODMCParticle*)arrayMC->At(ip);
 678     if(TMath::Abs(part->GetPdgCode())==thePDG){
 679       Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,part,fGoUpToQuark);
 680       fHistCheckOrigin->Fill(orig);
 681       if(fPromptFeeddown==kFeeddown && orig!=5) continue;
 682       else if(fPromptFeeddown==kPrompt && orig!=4) continue;
 683       else if(fPromptFeeddown==kBoth && orig<4) continue;
 684       fHistCheckOriginSel->Fill(orig);
 685       Int_t deca=0;
 686       Bool_t isGoodDecay=kFALSE;
 687       Int_t labDau[4]={-1,-1,-1,-1};
 688       if(fMeson==kDzero){
 689         deca=AliVertexingHFUtils::CheckD0Decay(arrayMC,part,labDau);
 690         if(part->GetNDaughters()!=2) continue;
 691         if(deca==1) isGoodDecay=kTRUE;
 692       }else if(fMeson==kDplus){
 693         deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,part,labDau);
 694         if(deca>0) isGoodDecay=kTRUE;
 695       }
 696       fHistCheckDecChan->Fill(deca);
 697       if(labDau[0]==-1){
 698         //      printf(Form("Meson %d Label of daughters not filled correctly -- %d\n",fMeson,isGoodDecay));
 699         continue; //protection against unfilled array of labels
 700       }
 701       Bool_t isInAcc=CheckAcceptance(arrayMC,nProng,labDau);
 702       if(isInAcc) fHistCheckDecChanAcc->Fill(deca);
 703       if(isGoodDecay){
 704         Double_t ptgen=part->Pt();
 705         Double_t ygen=part->Y();
 706         if(fAnalysisCuts->IsInFiducialAcceptance(ptgen,ygen)){
 707           fPtVsYGen->Fill(ptgen,ygen);
 708           if(TMath::Abs(ygen)<0.5) fPtVsYGenLimAcc->Fill(ptgen,ygen);
 709           if(isInAcc) fPtVsYGenAcc->Fill(ptgen,ygen);
 710         }
 711         if(TMath::Abs(ygen)<0.9) fPtVsYGenLargeAcc->Fill(ptgen,ygen);
 712       }
 713     }
 714   }
 715 }
 716
 717 //________________________________________________________________________
 718 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){
 719   // Fill histos for candidates with proper charge sign
 720
 721   Bool_t accept=kFALSE;
 722
 723   tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
 724   Double_t pt = tmpRD->Pt();
 725   Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
 726   Double_t mass=TMath::Sqrt(minv2);
 727
 728   if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
 729     Double_t rapid = tmpRD->Y(pdgD);
 730     if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
 731       fMassVsPtVsY->Fill(mass,pt,rapid);
 732       accept=kTRUE;
 733       if(fReadMC){
 734         Int_t signPdg[3]={0,0,0};
 735         for(Int_t iii=0; iii<nProngs; iii++) signPdg[iii]=pdgdau[iii];
 736         Int_t labD = tmpRD->MatchToMC(pdgD,arrayMC,nProngs,signPdg);
 737         if(labD>=0){
 738           AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(TMath::Abs(dgLabels[0])));
 739           if(part){
 740             Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,part,fGoUpToQuark);
 741             if((fPromptFeeddown==kFeeddown && orig==5)|| (fPromptFeeddown==kPrompt && orig==4) || (fPromptFeeddown==kBoth && orig>=4)) {
 742
 743               Int_t pdgCode = TMath::Abs( part->GetPdgCode() );
 744               if(pdgCode==321){
 745                 fMassVsPtVsYSig->Fill(mass,pt,rapid);
 746                 AliAODMCParticle* dmes =  dynamic_cast<AliAODMCParticle*>(arrayMC->At(labD));
 747                 fPtVsYReco->Fill(dmes->Pt(),dmes->Y());
 748               }else{
 749                 fMassVsPtVsYRefl->Fill(mass,pt,rapid);
 750               }
 751             }
 752           }
 753         }else{
 754           fMassVsPtVsYBkg->Fill(mass,pt,rapid);
 755         }
 756       }
 757     }
 758   }
 759
 760   Int_t nRotated=0;
 761   Double_t massRot=0;// calculated later only if candidate is acceptable
 762   Double_t angleProngXY;
 763   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])));
 764   else {
 765     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])));
 766   }
 767   Double_t ptOrig=pt;
 768
 769
 770   Double_t rotStep=(fMaxAngleForRot-fMinAngleForRot)/(fNRotations-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot
 771   Double_t rotStep3=(fMaxAngleForRot3-fMinAngleForRot3)/(fNRotations-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot
 772
 773   for(Int_t irot=0; irot<fNRotations; irot++){
 774     Double_t phirot=fMinAngleForRot+rotStep*irot;
 775     Double_t tmpx=px[0];
 776     Double_t tmpy=py[0];
 777     Double_t tmpx2=px[2];
 778     Double_t tmpy2=py[2];
 779     px[0]=tmpx*TMath::Cos(phirot)-tmpy*TMath::Sin(phirot);
 780     py[0]=tmpx*TMath::Sin(phirot)+tmpy*TMath::Cos(phirot);
 781     if(pdgD==411 || pdgD==431){
 782       Double_t phirot2=fMinAngleForRot3+rotStep3*irot;
 783       px[2]=tmpx*TMath::Cos(phirot2)-tmpy*TMath::Sin(phirot2);
 784       py[2]=tmpx*TMath::Sin(phirot2)+tmpy*TMath::Cos(phirot2);
 785     }
 786     tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
 787     pt = tmpRD->Pt();
 788     minv2 = tmpRD->InvMass2(nProngs,pdgdau);
 789     if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
 790       Double_t rapid = tmpRD->Y(pdgD);
 791       if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
 792         massRot=TMath::Sqrt(minv2);
 793         fMassVsPtVsYRot->Fill(massRot,pt,rapid);
 794         nRotated++;
 795         fDeltaMass->Fill(massRot-mass);
 796         if(fFullAnalysis){
 797           Double_t pointRot[5]={mass,massRot-mass,ptOrig,pt-ptOrig,angleProngXY};
 798           fDeltaMassFullAnalysis->Fill(pointRot);
 799         }
 800       }
 801     }
 802     px[0]=tmpx;
 803     py[0]=tmpy;
 804     if(pdgD==411 || pdgD==431){
 805       px[2]=tmpx2;
 806       py[2]=tmpy2;
 807     }
 808   }
 809   fNormRotated->Fill(nRotated);
 810
 811   return accept;
 812
 813 }
 814 //________________________________________________________________________
 815 void AliAnalysisTaskCombinHF::FillMEHistos(Int_t pdgD,Int_t nProngs, AliAODRecoDecay* tmpRD, Double_t* px, Double_t* py, Double_t* pz, UInt_t *pdgdau){
 816   // Fill histos for candidates in MixedEvents
 817
 818   tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
 819   Double_t pt = tmpRD->Pt();
 820   Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
 821   Double_t mass=TMath::Sqrt(minv2);
 822
 823   if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
 824     Double_t rapid = tmpRD->Y(pdgD);
 825     if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
 826       fMassVsPtVsYME->Fill(mass,pt,rapid);
 827     }
 828   }
 829   return;
 830 }
 831
 832 //________________________________________________________________________
 833 Bool_t AliAnalysisTaskCombinHF::IsTrackSelected(AliAODTrack* track){
 834   // track selection cuts
 835
 836   if(track->Charge()==0) return kFALSE;
 837   if(track->GetID()<0&&!fKeepNegID)return kFALSE;
 838   if(!(track->TestFilterMask(fFilterMask))) return kFALSE;
 839   if(!SelectAODTrack(track,fTrackCutsAll)) return kFALSE;
 840   return kTRUE;
 841 }
 842
 843 //________________________________________________________________________
 844 Bool_t AliAnalysisTaskCombinHF::IsKaon(AliAODTrack* track){
 845   // kaon selection cuts
 846
 847   if(!fPidHF) return kTRUE;
 848   Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);
 849   Double_t mom=track->P();
 850   if(SelectAODTrack(track,fTrackCutsKaon)) {
 851     if(isKaon>=1)    return kTRUE;
 852     if(isKaon<=-1) return kFALSE;
 853     switch(fPIDselCaseZero){// isKaon=0
 854       case 0:
 855       {
 856         return kTRUE;// always accept
 857       }
 858         break;
 859       case 1:
 860       {
 861         if(isKaon>=0 && track->P()>fmaxPforIDKaon)return kTRUE;// accept only if in a compatibility band starting from p=fmaxPforIDKaon
 862       }
 863         break;
 864       case 2:
 865       {
 866         if(track->P()>fmaxPforIDKaon)return kTRUE;
 867         AliPIDResponse *pidResp=fPidHF->GetPidResponse();// more elaborated strategy: asymmetric cuts, with fix momenta and nsigma ranges for the moment
 868         Double_t nsigma=pidResp->NumberOfSigmasTPC(track,AliPID::kKaon);
 869         if(nsigma>-2.&& nsigma<3. && mom<0.6)isKaon=1;
 870         else if(nsigma>-1.&& nsigma<3.&& mom<0.8)isKaon=1;
 871         if(isKaon==1)return kTRUE;
 872       }
 873         break;
 874       default:
 875       {
 876         AliWarning(Form("WRONG CASE OF PID STRATEGY SELECTED: %d (can range from 0 to 2)",fPIDselCaseZero));
 877         return kFALSE;// actually case 0 could be set as the default and return kTRUE
 878       }
 879     }
 880   }
 881
 882   return kFALSE;
 883 }
 884 //_______________________________________________________________________
 885 Bool_t AliAnalysisTaskCombinHF::IsPion(AliAODTrack* track){
 886   // pion selection cuts
 887
 888   if(!fPidHF) return kTRUE;
 889   Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
 890   Double_t mom=track->P();
 891   if(SelectAODTrack(track,fTrackCutsPion)) {
 892     if(isPion>=1)    return kTRUE;
 893     if(isPion<=-1) return kFALSE;
 894     switch(fPIDselCaseZero){// isPion=0
 895       case 0:
 896       {
 897         return kTRUE;// always accept
 898       }
 899         break;
 900       case 1:
 901       {
 902         if(track->P()>fmaxPforIDPion)return kTRUE;// accept only if in a compatibility band starting from p=fmaxPforIDPion
 903       }
 904         break;
 905       case 2:
 906       {
 907         // more elaborated strategy: asymmetric cuts, with fix momenta and nsigma ranges for the moment
 908         if(track->P()>fmaxPforIDPion)return kTRUE;
 909         AliPIDResponse *pidResp=fPidHF->GetPidResponse();
 910         Double_t nsigma=pidResp->NumberOfSigmasTPC(track,AliPID::kPion);
 911         if(nsigma<2.&& nsigma>-3. && mom<0.6)isPion=1;
 912         else if(nsigma<1. && nsigma> -3. && mom<0.8)isPion=1;
 913         if(isPion==1)return kTRUE;
 914       }
 915         break;
 916       default:
 917       {
 918         AliWarning(Form("WRONG CASE OF PID STRATEGY SELECTED: %d (can range from 0 to 2)",fPIDselCaseZero));
 919         return kFALSE;// actually case 0 could be set as the default and return kTRUE
 920       }
 921     }
 922   }
 923
 924   return kFALSE;
 925 }
 926
 927 //________________________________________________________________________
 928 Bool_t AliAnalysisTaskCombinHF::SelectAODTrack(AliAODTrack *track, AliESDtrackCuts *cuts){
 929   // AOD track selection
 930
 931   if(!cuts) return kTRUE;
 932
 933   AliESDtrack esdTrack(track);
 934   // set the TPC cluster info
 935   esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
 936   esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
 937   esdTrack.SetTPCPointsF(track->GetTPCNclsF());
 938   if(!cuts->IsSelected(&esdTrack)) return kFALSE;
 939
 940   return kTRUE;
 941 }
 942
 943 //_________________________________________________________________
 944 Bool_t AliAnalysisTaskCombinHF::CheckAcceptance(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
 945   // check if the decay products are in the good eta and pt range
 946   for (Int_t iProng = 0; iProng<nProng; iProng++){
 947     AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
 948     if(!mcPartDaughter) return kFALSE;
 949     Double_t eta = mcPartDaughter->Eta();
 950     Double_t pt = mcPartDaughter->Pt();
 951     if (TMath::Abs(eta) > fEtaAccCut || pt < fPtAccCut) return kFALSE;
 952   }
 953   return kTRUE;
 954 }
 955 //_________________________________________________________________
 956 Bool_t AliAnalysisTaskCombinHF::CanBeMixed(Double_t zv1, Double_t zv2, Double_t mult1, Double_t mult2){
 957   if(fUsePoolsZ && fzVertPoolLims){
 958     Int_t theBin1=TMath::BinarySearch(fNzVertPoolsLimSize,fzVertPoolLims,zv1);
 959     Int_t theBin2=TMath::BinarySearch(fNzVertPoolsLimSize,fzVertPoolLims,zv2);
 960     if(theBin1!=theBin2) return kFALSE;
 961     if(theBin1<0 || theBin2<0) return kFALSE;
 962     if(theBin1>=fNzVertPools || theBin2>=fNzVertPools) return kFALSE;
 963   }else{
 964     if(TMath::Abs(zv2-zv1)>fMaxzVertDistForMix) return kFALSE;
 965   }
 966
 967   if(fUsePoolsM && fMultPoolLims){
 968     Int_t theBin1=TMath::BinarySearch(fNMultPoolsLimSize,fMultPoolLims,mult1);
 969     Int_t theBin2=TMath::BinarySearch(fNMultPoolsLimSize,fMultPoolLims,mult2);
 970     if(theBin1!=theBin2) return kFALSE;
 971     if(theBin1<0 || theBin2<0) return kFALSE;
 972     if(theBin1>=fNMultPools || theBin2>=fNMultPools) return kFALSE;
 973   }else{
 974     if(TMath::Abs(mult2-mult1)>fMaxMultDiffForMix) return kFALSE;
 975   }
 976   return kTRUE;
 977 }
 978
 979 //_________________________________________________________________
 980 void AliAnalysisTaskCombinHF::FinishTaskOutput()
 981 {
 982   // perform mixed event analysis
 983   if(!fDoEventMixing) return;
 984
 985   Int_t nEvents=fEventBuffer->GetEntries();
 986   printf("Start Event Mixing of %d events\n",nEvents);
 987
 988   TObjArray* karray=0x0;
 989   TObjArray* parray=0x0;
 990   Double_t zVertex,mult;
 991   fEventBuffer->SetBranchAddress("karray", &karray);
 992   fEventBuffer->SetBranchAddress("parray", &parray);
 993   fEventBuffer->SetBranchAddress("zVertex", &zVertex);
 994   fEventBuffer->SetBranchAddress("multiplicity", &mult);
 995
 996   // dummy values of track impact parameter, needed to build an AliAODRecoDecay object
 997   Double_t d02[2]={0.,0.};
 998   Double_t d03[3]={0.,0.,0.};
 999   AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
1000   AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
1001   UInt_t pdg0[2]={321,211};
1002   UInt_t pdgp[3]={321,211,211};
1003   Double_t px[3],py[3],pz[3];
1004
1005   for(Int_t iEv1=0; iEv1<nEvents; iEv1++){
1006     fEventBuffer->GetEvent(iEv1);
1007     Double_t zVertex1=zVertex;
1008     Double_t mult1=mult;
1009     TObjArray* karray1=new TObjArray(*karray);
1010     for(Int_t iEv2=0; iEv2<fMaxNumberOfEventsForMixing; iEv2++){
1011       Int_t iToMix=iEv1+iEv2+1;
1012       if(iEv1>=(nEvents-fMaxNumberOfEventsForMixing)) iToMix=iEv1-iEv2-1;
1013       if(iToMix<0) continue;
1014       if(iToMix==iEv1) continue;
1015       fEventBuffer->GetEvent(iToMix);
1016       Double_t zVertex2=zVertex;
1017       Double_t mult2=mult;
1018       if(CanBeMixed(zVertex1,zVertex2,mult1,mult2)){
1019         TObjArray* parray2=new TObjArray(*parray);
1020         Int_t nKaons=karray1->GetEntries();
1021         Int_t nPions=parray2->GetEntries();
1022         Bool_t doThird=kTRUE;
1023         Int_t iToMix3=iEv1+2*fMaxNumberOfEventsForMixing+1;
1024         if(iToMix3>=nEvents) iToMix3=iEv1-2*fMaxNumberOfEventsForMixing-1;
1025         if(iToMix3<0 || iToMix3==iEv1 || iToMix3==iToMix) doThird=kFALSE;
1026         TObjArray* parray3=0x0;
1027         Double_t zVertex3=-999.;
1028         Double_t mult3=999999;
1029         if(fMeson!=kDzero && doThird){
1030           fEventBuffer->GetEvent(iToMix3);
1031           zVertex3=zVertex;
1032           mult3=mult;
1033           if(CanBeMixed(zVertex1,zVertex3,mult1,mult3)){
1034             parray3=new TObjArray(*parray);
1035           }else{
1036             doThird=kFALSE;
1037           }
1038         }
1039
1040         for(Int_t iTr1=0; iTr1<nKaons; iTr1++){
1041           TLorentzVector* trK=(TLorentzVector*)karray1->At(iTr1);
1042           Double_t chargeK=trK->T();
1043           px[0] = trK->Px();
1044           py[0] = trK->Py();
1045           pz[0] = trK->Pz();
1046           for(Int_t iTr2=0; iTr2<nPions; iTr2++){
1047             TLorentzVector* trPi1=(TLorentzVector*)parray2->At(iTr2);
1048             Double_t chargePi1=trPi1->T();
1049             px[1] = trPi1->Px();
1050             py[1] = trPi1->Py();
1051             pz[1] = trPi1->Pz();
1052             if(chargePi1*chargeK<0){
1053               if(fMeson==kDzero){
1054                 FillMEHistos(421,2,tmpRD2,px,py,pz,pdg0);
1055               }else{
1056                 if(doThird && parray3){
1057                   Int_t nPions3=parray3->GetEntries();
1058                   for(Int_t iTr3=iTr2+1; iTr3<nPions3; iTr3++){
1059                     TLorentzVector* trPi2=(TLorentzVector*)parray3->At(iTr3);
1060                     Double_t chargePi2=trPi2->T();
1061                     px[2] = trPi2->Px();
1062                     py[2] = trPi2->Py();
1063                     pz[2] = trPi2->Pz();
1064                     if(chargePi2*chargeK<0){
1065                       FillMEHistos(411,3,tmpRD3,px,py,pz,pdgp);
1066                     }
1067                   }
1068                 }
1069               }
1070             }
1071           }
1072         }
1073         delete parray3;
1074         delete parray2;
1075       }
1076     }
1077     delete karray1;
1078   }
1079   delete tmpRD2;
1080   delete tmpRD3;
1081 }
1082 //_________________________________________________________________
1083 void AliAnalysisTaskCombinHF::Terminate(Option_t */*option*/)
1084 {
1085   // Terminate analysis
1086   //
1087   if(fDebug > 1) printf("AliAnalysisTaskCombinHF: Terminate() \n");
1088   fOutput = dynamic_cast<TList*> (GetOutputData(1));
1089   if (!fOutput) {
1090     printf("ERROR: fOutput not available\n");
1091     return;
1092   }
1093   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
1094   if(fHistNEvents){
1095     printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
1096   }else{
1097     printf("ERROR: fHistNEvents not available\n");
1098     return;
1099   }
1100   return;
1101 }
1102