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