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   fFilterMask(BIT(4)),
  73   fTrackCutsAll(0x0),
  74   fTrackCutsPion(0x0),
  75   fTrackCutsKaon(0x0),
  76   fPidHF(new AliAODPidHF()),
  77   fAnalysisCuts(0x0),
  78   fMinMass(1.720),
  79   fMaxMass(2.150),
  80   fMaxPt(10.),
  81   fEtaAccCut(0.9),
  82   fPtAccCut(0.1),
  83   fNRotations(9),
  84   fMinAngleForRot(5*TMath::Pi()/6),
  85   fMaxAngleForRot(7*TMath::Pi()/6),
  86   fMinAngleForRot3(2*TMath::Pi()/6),
  87   fMaxAngleForRot3(4*TMath::Pi()/6),
  88   fCounter(0x0),
  89   fMeson(kDzero),
  90   fReadMC(kFALSE),
  91   fPromptFeeddown(kPrompt),
  92   fGoUpToQuark(kTRUE),
  93   fFullAnalysis(0),
  94   fPIDstrategy(knSigma),
  95   fmaxPforIDPion(0.8),
  96   fmaxPforIDKaon(2.),
  97   fKeepNegID(kFALSE),
  98   fPIDselCaseZero(0),
  99   fBayesThresKaon(0.4),
 100   fBayesThresPion(0.4)
 101 {
 102   // default constructor
 103 }
 104
 105 //________________________________________________________________________
 106 AliAnalysisTaskCombinHF::AliAnalysisTaskCombinHF(Int_t meson, AliRDHFCuts* analysiscuts):
 107   AliAnalysisTaskSE("DmesonCombin"),
 108   fOutput(0x0),
 109   fHistNEvents(0x0),
 110   fHistTrackStatus(0x0),
 111   fHistCheckOrigin(0x0),
 112   fHistCheckOriginSel(0x0),
 113   fHistCheckDecChan(0x0),
 114   fHistCheckDecChanAcc(0x0),
 115   fPtVsYGen(0x0),
 116   fPtVsYGenLargeAcc(0x0),
 117   fPtVsYGenLimAcc(0x0),
 118   fPtVsYGenAcc(0x0),
 119   fPtVsYReco(0x0),
 120   fMassVsPtVsY(0x0),
 121   fMassVsPtVsYRot(0x0),
 122   fMassVsPtVsYLSpp(0x0),
 123   fMassVsPtVsYLSmm(0x0),
 124   fMassVsPtVsYSig(0x0),
 125   fMassVsPtVsYRefl(0x0),
 126   fMassVsPtVsYBkg(0x0),
 127   fNSelected(0x0),
 128   fNormRotated(0x0),
 129   fDeltaMass(0x0),
 130   fDeltaMassFullAnalysis(0x0),
 131   fFilterMask(BIT(4)),
 132   fTrackCutsAll(0x0),
 133   fTrackCutsPion(0x0),
 134   fTrackCutsKaon(0x0),
 135   fPidHF(new AliAODPidHF()),
 136   fAnalysisCuts(analysiscuts),
 137   fMinMass(1.720),
 138   fMaxMass(2.150),
 139   fMaxPt(10.),
 140   fEtaAccCut(0.9),
 141   fPtAccCut(0.1),
 142   fNRotations(9),
 143   fMinAngleForRot(5*TMath::Pi()/6),
 144   fMaxAngleForRot(7*TMath::Pi()/6),
 145   fMinAngleForRot3(2*TMath::Pi()/6),
 146   fMaxAngleForRot3(4*TMath::Pi()/6),
 147   fCounter(0x0),
 148   fMeson(meson),
 149   fReadMC(kFALSE),
 150   fPromptFeeddown(1),
 151   fGoUpToQuark(kTRUE),
 152   fFullAnalysis(0),
 153   fPIDstrategy(knSigma),
 154   fmaxPforIDPion(0.8),
 155   fmaxPforIDKaon(2.),
 156   fKeepNegID(kFALSE),
 157   fPIDselCaseZero(0),
 158   fBayesThresKaon(0.4),
 159   fBayesThresPion(0.4)
 160 {
 161   // standard constructor
 162
 163   DefineOutput(1,TList::Class());  //My private output
 164   DefineOutput(2,AliNormalizationCounter::Class());
 165 }
 166
 167 //________________________________________________________________________
 168 AliAnalysisTaskCombinHF::~AliAnalysisTaskCombinHF()
 169 {
 170   //
 171   // Destructor
 172   //
 173   delete fOutput;
 174   delete fHistNEvents;
 175   delete fHistTrackStatus;
 176   delete fHistCheckOrigin;
 177   delete fHistCheckOriginSel;
 178   delete fHistCheckDecChan;
 179   delete fHistCheckDecChanAcc;
 180   delete fPtVsYGen;
 181   delete fPtVsYGenLargeAcc;
 182   delete fPtVsYGenLimAcc;
 183   delete fPtVsYGenAcc;
 184   delete fPtVsYReco;
 185   delete fMassVsPtVsY;
 186   delete fMassVsPtVsYLSpp;
 187   delete fMassVsPtVsYLSmm;
 188   delete fMassVsPtVsYRot;
 189   delete fMassVsPtVsYSig;
 190   delete fMassVsPtVsYRefl;
 191   delete fMassVsPtVsYBkg;
 192   delete fNSelected;
 193   delete fNormRotated;
 194   delete fDeltaMass;
 195   delete fCounter;
 196   delete fTrackCutsAll;
 197   delete fTrackCutsPion;
 198   delete fTrackCutsKaon;
 199   delete fPidHF;
 200   delete fAnalysisCuts;
 201 }
 202
 203 //________________________________________________________________________
 204 void AliAnalysisTaskCombinHF::UserCreateOutputObjects()
 205 {
 206   // Create the output container
 207   //
 208   if(fDebug > 1) printf("AnalysisTaskCombinHF::UserCreateOutputObjects() \n");
 209
 210   fOutput = new TList();
 211   fOutput->SetOwner();
 212   fOutput->SetName("OutputHistos");
 213
 214   fHistNEvents = new TH1F("hNEvents", "number of events ",8,-0.5,7.5);
 215   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 216   fHistNEvents->GetXaxis()->SetBinLabel(2,"n. passing IsEvSelected");
 217   fHistNEvents->GetXaxis()->SetBinLabel(3,"n. rejected due to trigger");
 218   fHistNEvents->GetXaxis()->SetBinLabel(4,"n. rejected due to not reco vertex");
 219   fHistNEvents->GetXaxis()->SetBinLabel(5,"n. rejected for contr vertex");
 220   fHistNEvents->GetXaxis()->SetBinLabel(6,"n. rejected for vertex out of accept");
 221   fHistNEvents->GetXaxis()->SetBinLabel(7,"n. rejected for pileup events");
 222   fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events");
 223
 224   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
 225   fHistNEvents->Sumw2();
 226   fHistNEvents->SetMinimum(0);
 227   fOutput->Add(fHistNEvents);
 228
 229   fHistTrackStatus  = new TH1F("hTrackStatus", "",8,-0.5,7.5);
 230   fHistTrackStatus->GetXaxis()->SetBinLabel(1,"Not OK");
 231   fHistTrackStatus->GetXaxis()->SetBinLabel(2,"Track OK");
 232   fHistTrackStatus->GetXaxis()->SetBinLabel(3,"Kaon, Not OK");
 233   fHistTrackStatus->GetXaxis()->SetBinLabel(4,"Kaon OK");
 234   fHistTrackStatus->GetXaxis()->SetBinLabel(5,"Pion, Not OK");
 235   fHistTrackStatus->GetXaxis()->SetBinLabel(6,"Pion OK");
 236   fHistTrackStatus->GetXaxis()->SetBinLabel(7,"Kaon||Pion, Not OK");
 237   fHistTrackStatus->GetXaxis()->SetBinLabel(8,"Kaon||Pion OK");
 238   fHistTrackStatus->GetXaxis()->SetNdivisions(1,kFALSE);
 239   fHistTrackStatus->Sumw2();
 240   fHistTrackStatus->SetMinimum(0);
 241   fOutput->Add(fHistTrackStatus);
 242
 243   Int_t nPtBins = fMaxPt*10;
 244
 245   if(fReadMC){
 246
 247     fHistCheckOrigin=new TH1F("hCheckOrigin","",7,-1.5,5.5);
 248     fHistCheckOrigin->Sumw2();
 249     fHistCheckOrigin->SetMinimum(0);
 250     fOutput->Add(fHistCheckOrigin);
 251
 252     fHistCheckOriginSel=new TH1F("hCheckOriginSel","",7,-1.5,5.5);
 253     fHistCheckOriginSel->Sumw2();
 254     fHistCheckOriginSel->SetMinimum(0);
 255     fOutput->Add(fHistCheckOriginSel);
 256
 257     fHistCheckDecChan=new TH1F("hCheckDecChan","",7,-2.5,4.5);
 258     fHistCheckDecChan->Sumw2();
 259     fHistCheckDecChan->SetMinimum(0);
 260     fOutput->Add(fHistCheckDecChan);
 261
 262     fHistCheckDecChanAcc=new TH1F("hCheckDecChanAcc","",7,-2.5,4.5);
 263     fHistCheckDecChanAcc->Sumw2();
 264     fHistCheckDecChanAcc->SetMinimum(0);
 265     fOutput->Add(fHistCheckDecChanAcc);
 266
 267     fPtVsYGen= new TH2F("hPtVsYGen","",nPtBins,0.,fMaxPt,20,-1.,1.);
 268     fPtVsYGen->Sumw2();
 269     fPtVsYGen->SetMinimum(0);
 270     fOutput->Add(fPtVsYGen);
 271
 272     fPtVsYGenLargeAcc= new TH2F("hPtVsYGenLargeAcc","",nPtBins,0.,fMaxPt,20,-1.,1.);
 273     fPtVsYGenLargeAcc->Sumw2();
 274     fPtVsYGenLargeAcc->SetMinimum(0);
 275     fOutput->Add(fPtVsYGenLargeAcc);
 276
 277     fPtVsYGenLimAcc= new TH2F("hPtVsYGenLimAcc","",nPtBins,0.,fMaxPt,20,-1.,1.);
 278     fPtVsYGenLimAcc->Sumw2();
 279     fPtVsYGenLimAcc->SetMinimum(0);
 280     fOutput->Add(fPtVsYGenLimAcc);
 281
 282     fPtVsYGenAcc= new TH2F("hPtVsYGenAcc","",nPtBins,0.,fMaxPt,20,-1.,1.);
 283     fPtVsYGenAcc->Sumw2();
 284     fPtVsYGenAcc->SetMinimum(0);
 285     fOutput->Add(fPtVsYGenAcc);
 286
 287     fPtVsYReco= new TH2F("hPtVsYReco","",nPtBins,0.,fMaxPt,20,-1.,1.);
 288     fPtVsYReco->Sumw2();
 289     fPtVsYReco->SetMinimum(0);
 290     fOutput->Add(fPtVsYReco);
 291   }
 292
 293
 294   Int_t nMassBins=fMaxMass*1000.-fMinMass*1000.;
 295   Double_t maxm=fMinMass+nMassBins*0.001;
 296   fMassVsPtVsY=new TH3F("hMassVsPtVsY","",nMassBins,fMinMass,maxm,nPtBins,0.,fMaxPt,20,-1.,1.);
 297   fMassVsPtVsY->Sumw2();
 298   fMassVsPtVsY->SetMinimum(0);
 299   fOutput->Add(fMassVsPtVsY);
 300
 301   fMassVsPtVsYRot=new TH3F("hMassVsPtVsYRot","",nMassBins,fMinMass,maxm,nPtBins,0.,fMaxPt,20,-1.,1.);
 302   fMassVsPtVsYRot->Sumw2();
 303   fMassVsPtVsYRot->SetMinimum(0);
 304   fOutput->Add(fMassVsPtVsYRot);
 305
 306   if(fMeson==kDzero){
 307     fMassVsPtVsYLSpp=new TH3F("hMassVsPtVsYLSpp","",nMassBins,fMinMass,maxm,nPtBins,0.,fMaxPt,20,-1.,1.);
 308     fMassVsPtVsYLSpp->Sumw2();
 309     fMassVsPtVsYLSpp->SetMinimum(0);
 310     fOutput->Add(fMassVsPtVsYLSpp);
 311     fMassVsPtVsYLSmm=new TH3F("hMassVsPtVsYLSmm","",nMassBins,fMinMass,maxm,nPtBins,0.,fMaxPt,20,-1.,1.);
 312     fMassVsPtVsYLSmm->Sumw2();
 313     fMassVsPtVsYLSmm->SetMinimum(0);
 314     fOutput->Add(fMassVsPtVsYLSmm);
 315   }
 316
 317   fMassVsPtVsYSig=new TH3F("hMassVsPtVsYSig","",nMassBins,fMinMass,maxm,nPtBins,0.,fMaxPt,20,-1.,1.);
 318   fMassVsPtVsYSig->Sumw2();
 319   fMassVsPtVsYSig->SetMinimum(0);
 320   fOutput->Add(fMassVsPtVsYSig);
 321
 322   fMassVsPtVsYRefl=new TH3F("hMassVsPtVsYRefl","",nMassBins,fMinMass,maxm,nPtBins,0.,fMaxPt,20,-1.,1.);
 323   fMassVsPtVsYRefl->Sumw2();
 324   fMassVsPtVsYRefl->SetMinimum(0);
 325   fOutput->Add(fMassVsPtVsYRefl);
 326
 327   fMassVsPtVsYBkg=new TH3F("hMassVsPtVsYBkg","",nMassBins,fMinMass,maxm,nPtBins,0.,fMaxPt,20,-1.,1.);
 328   fMassVsPtVsYBkg->Sumw2();
 329   fMassVsPtVsYBkg->SetMinimum(0);
 330   fOutput->Add(fMassVsPtVsYBkg);
 331
 332   fNSelected=new TH1F("hNSelected","",100,-0.5,99.5);
 333   fNSelected->Sumw2();
 334   fNSelected->SetMinimum(0);
 335   fOutput->Add(fNSelected);
 336
 337   fNormRotated=new TH1F("hNormRotated","",11,-0.5,10.5);
 338   fNormRotated->Sumw2();
 339   fNormRotated->SetMinimum(0);
 340   fOutput->Add(fNormRotated);
 341
 342   fDeltaMass=new TH1F("hDeltaMass","",100,-0.4,0.4);
 343   fDeltaMass->Sumw2();
 344   fDeltaMass->SetMinimum(0);
 345   fOutput->Add(fDeltaMass);
 346
 347   Int_t binSparseDMassRot[5]={nMassBins,100,24,40,20};
 348   Double_t edgeLowSparseDMassRot[5]={fMinMass,-0.4,0.,-4.,0};
 349   Double_t edgeHighSparseDMassRot[5]={maxm,0.4,12.,4.,3.14};
 350   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);
 351   fOutput->Add(fDeltaMassFullAnalysis);
 352
 353   //Counter for Normalization
 354   fCounter = new AliNormalizationCounter("NormalizationCounter");
 355   fCounter->Init();
 356
 357   PostData(1,fOutput);
 358   PostData(2,fCounter);
 359 }
 360
 361 //________________________________________________________________________
 362 void AliAnalysisTaskCombinHF::UserExec(Option_t */*option*/){
 363   //Build the 3-track combinatorics (+-+ and -+-) for D+->Kpipi decays
 364
 365   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 366   if(!aod && AODEvent() && IsStandardAOD()) {
 367     // In case there is an AOD handler writing a standard AOD, use the AOD
 368     // event in memory rather than the input (ESD) event.
 369     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 370   }
 371   if(!aod){
 372     printf("AliAnalysisTaskCombinHF::UserExec: AOD not found!\n");
 373     return;
 374   }
 375
 376   // fix for temporary bug in ESDfilter
 377   // the AODs with null vertex pointer didn't pass the PhysSel
 378   if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
 379   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 380   AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
 381   AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
 382   fPidHF->SetPidResponse(pidResp);
 383
 384
 385   fHistNEvents->Fill(0); // count event
 386   // Post the data already here
 387   PostData(1,fOutput);
 388
 389   fCounter->StoreEvent(aod,fAnalysisCuts,fReadMC);
 390
 391   Bool_t isEvSel=fAnalysisCuts->IsEventSelected(aod);
 392   if(fAnalysisCuts->IsEventRejectedDueToTrigger())fHistNEvents->Fill(2);
 393   if(fAnalysisCuts->IsEventRejectedDueToNotRecoVertex())fHistNEvents->Fill(3);
 394   if(fAnalysisCuts->IsEventRejectedDueToVertexContributors())fHistNEvents->Fill(4);
 395   if(fAnalysisCuts->IsEventRejectedDueToZVertexOutsideFiducialRegion())fHistNEvents->Fill(5);
 396   //  if(fAnalysisCuts->IsEventRejectedDueToPileup())fHistNEvents->Fill(6);
 397   if(fAnalysisCuts->IsEventRejectedDueToCentrality()) fHistNEvents->Fill(7);
 398
 399   if(!isEvSel)return;
 400
 401   fHistNEvents->Fill(1);
 402
 403   TClonesArray *arrayMC=0;
 404   AliAODMCHeader *mcHeader=0;
 405   if(fReadMC){
 406     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 407     if(!arrayMC) {
 408       printf("AliAnalysisTaskCombinHF::UserExec: MC particles branch not found!\n");
 409       return;
 410     }
 411
 412     // load MC header
 413     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 414     if(!mcHeader) {
 415       printf("AliAnalysisTaskCombinHF::UserExec: MC header branch not found!\n");
 416       return;
 417     }
 418     FillGenHistos(arrayMC);
 419   }
 420
 421
 422   Int_t ntracks=aod->GetNTracks();
 423
 424   // select and flag tracks
 425   UChar_t* status = new UChar_t[ntracks];
 426   for(Int_t iTr=0; iTr<ntracks; iTr++){
 427     status[iTr]=0;
 428     AliAODTrack* track=aod->GetTrack(iTr);
 429     if(IsTrackSelected(track)) status[iTr]+=1;
 430
 431     // PID
 432     if (fPIDstrategy == knSigma) {
 433       // nsigma PID
 434       if(IsKaon(track)) status[iTr]+=2;
 435       if(IsPion(track)) status[iTr]+=4;
 436     }
 437     else if (fPIDstrategy == kBayesianMaxProb || fPIDstrategy == kBayesianThres) {
 438       // Bayesian PID
 439       Double_t *weights = new Double_t[AliPID::kSPECIES];
 440       fPidHF->GetPidCombined()->ComputeProbabilities(track, fPidHF->GetPidResponse(), weights);
 441       if (fPIDstrategy == kBayesianMaxProb) {
 442         if (TMath::MaxElement(AliPID::kSPECIES, weights) == weights[AliPID::kKaon]) status[iTr] += 2;
 443         if (TMath::MaxElement(AliPID::kSPECIES, weights) == weights[AliPID::kPion]) status[iTr] += 4;
 444       }
 445       if (fPIDstrategy == kBayesianThres) {
 446         if (weights[AliPID::kKaon] > fBayesThresKaon) status[iTr] += 2;
 447         if (weights[AliPID::kPion] > fBayesThresPion) status[iTr] += 4;
 448       }
 449       delete[] weights;
 450     }
 451
 452     fHistTrackStatus->Fill(status[iTr]);
 453   }
 454
 455   // build the combinatorics
 456   Int_t nSelected=0;
 457   Int_t nFiltered=0;
 458   Double_t dummypos[3]={0.,0.,0.};
 459   AliAODVertex* v2=new AliAODVertex(dummypos,999.,-1,2);
 460   AliAODVertex* v3=new AliAODVertex(dummypos,999.,-1,3);
 461   // dummy values of track impact parameter, needed to build an AliAODRecoDecay object
 462   Double_t d02[2]={0.,0.};
 463   Double_t d03[3]={0.,0.,0.};
 464   AliAODRecoDecay* tmpRD2 = new AliAODRecoDecay(0x0,2,0,d02);
 465   AliAODRecoDecay* tmpRD3 = new AliAODRecoDecay(0x0,3,1,d03);
 466   UInt_t pdg0[2]={321,211};
 467   UInt_t pdgp[3]={321,211,211};
 468   //  UInt_t pdgs[3]={321,321,211};
 469   Double_t tmpp[3];
 470   Double_t px[3],py[3],pz[3];
 471   Int_t dgLabels[3];
 472
 473   for(Int_t iTr1=0; iTr1<ntracks; iTr1++){
 474     AliAODTrack* trK=aod->GetTrack(iTr1);
 475     if((status[iTr1] & 1)==0) continue;
 476     if((status[iTr1] & 2)==0) continue;
 477     Int_t chargeK=trK->Charge();
 478     trK->GetPxPyPz(tmpp);
 479     px[0] = tmpp[0];
 480     py[0] = tmpp[1];
 481     pz[0] = tmpp[2];
 482     dgLabels[0]=trK->GetLabel();
 483     for(Int_t iTr2=0; iTr2<ntracks; iTr2++){
 484       if((status[iTr2] & 1)==0) continue;
 485       if((status[iTr2] & 4)==0) continue;
 486       if(iTr1==iTr2) continue;
 487       AliAODTrack* trPi1=aod->GetTrack(iTr2);
 488       Int_t chargePi1=trPi1->Charge();
 489       trPi1->GetPxPyPz(tmpp);
 490       px[1] = tmpp[0];
 491       py[1] = tmpp[1];
 492       pz[1] = tmpp[2];
 493       dgLabels[1]=trPi1->GetLabel();
 494       if(chargePi1==chargeK){
 495         if(fMeson==kDzero) FillLSHistos(421,2,tmpRD2,px,py,pz,pdg0,chargePi1);
 496         continue;
 497       }
 498       if(fMeson==kDzero){
 499         nFiltered++;
 500         v2->AddDaughter(trK);
 501         v2->AddDaughter(trPi1);
 502         tmpRD2->SetSecondaryVtx(v2);
 503         Bool_t ok=FillHistos(421,2,tmpRD2,px,py,pz,pdg0,arrayMC,dgLabels);
 504         v2->RemoveDaughters();
 505         if(ok) nSelected++;
 506       }else{
 507         for(Int_t iTr3=iTr2+1; iTr3<ntracks; iTr3++){
 508           if((status[iTr3] & 1)==0) continue;
 509           if((status[iTr3] & 4)==0) continue;
 510           if(iTr1==iTr3) continue;
 511           AliAODTrack* trPi2=aod->GetTrack(iTr3);
 512           Int_t chargePi2=trPi2->Charge();
 513           if(chargePi2==chargeK) continue;
 514           nFiltered++;
 515           trPi2->GetPxPyPz(tmpp);
 516           px[2] = tmpp[0];
 517           py[2] = tmpp[1];
 518           pz[2] = tmpp[2];
 519           dgLabels[2]=trPi2->GetLabel();
 520           v3->AddDaughter(trK);
 521           v3->AddDaughter(trPi1);
 522           v3->AddDaughter(trPi2);
 523           tmpRD3->SetSecondaryVtx(v3);
 524           Bool_t ok=FillHistos(411,3,tmpRD3,px,py,pz,pdgp,arrayMC,dgLabels);
 525           v3->RemoveDaughters();
 526           if(ok) nSelected++;
 527         }
 528       }
 529     }
 530   }
 531
 532   delete [] status;
 533   delete v2;
 534   delete v3;
 535   delete tmpRD2;
 536   delete tmpRD3;
 537
 538   fNSelected->Fill(nSelected);
 539
 540   fCounter->StoreCandidates(aod,nFiltered,kTRUE);
 541   fCounter->StoreCandidates(aod,nSelected,kFALSE);
 542
 543   PostData(1,fOutput);
 544   PostData(2,fCounter);
 545
 546   return;
 547 }
 548
 549 //________________________________________________________________________
 550 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){
 551   // Fill histos for LS candidates
 552
 553   tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
 554   Double_t pt = tmpRD->Pt();
 555   Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
 556   if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
 557     Double_t rapid = tmpRD->Y(pdgD);
 558     if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
 559       if(charge>0) fMassVsPtVsYLSpp->Fill(TMath::Sqrt(minv2),pt,rapid);
 560       else fMassVsPtVsYLSmm->Fill(TMath::Sqrt(minv2),pt,rapid);
 561     }
 562   }
 563   return;
 564 }
 565
 566 //________________________________________________________________________
 567 void AliAnalysisTaskCombinHF::FillGenHistos(TClonesArray* arrayMC){
 568   // Fill histos with generated quantities
 569   Int_t totPart=arrayMC->GetEntriesFast();
 570   Int_t thePDG=411;
 571   Int_t nProng=3;
 572   if(fMeson==kDzero){
 573     thePDG=421;
 574     nProng=2;
 575   }
 576   for(Int_t ip=0; ip<totPart; ip++){
 577     AliAODMCParticle *part = (AliAODMCParticle*)arrayMC->At(ip);
 578     if(TMath::Abs(part->GetPdgCode())==thePDG){
 579       Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,part,fGoUpToQuark);
 580       fHistCheckOrigin->Fill(orig);
 581       if(fPromptFeeddown==kFeeddown && orig!=5) continue;
 582       else if(fPromptFeeddown==kPrompt && orig!=4) continue;
 583       else if(fPromptFeeddown==kBoth && orig<4) continue;
 584       fHistCheckOriginSel->Fill(orig);
 585       Int_t deca=0;
 586       Bool_t isGoodDecay=kFALSE;
 587       Int_t labDau[4]={-1,-1,-1,-1};
 588       if(fMeson==kDzero){
 589         deca=AliVertexingHFUtils::CheckD0Decay(arrayMC,part,labDau);
 590         if(part->GetNDaughters()!=2) continue;
 591         if(deca==1) isGoodDecay=kTRUE;
 592       }else if(fMeson==kDplus){
 593         deca=AliVertexingHFUtils::CheckDplusDecay(arrayMC,part,labDau);
 594         if(deca>0) isGoodDecay=kTRUE;
 595       }
 596       fHistCheckDecChan->Fill(deca);
 597       if(labDau[0]==-1){
 598         //      printf(Form("Meson %d Label of daughters not filled correctly -- %d\n",fMeson,isGoodDecay));
 599         continue; //protection against unfilled array of labels
 600       }
 601       Bool_t isInAcc=CheckAcceptance(arrayMC,nProng,labDau);
 602       if(isInAcc) fHistCheckDecChanAcc->Fill(deca);
 603       if(isGoodDecay){
 604         Double_t ptgen=part->Pt();
 605         Double_t ygen=part->Y();
 606         if(fAnalysisCuts->IsInFiducialAcceptance(ptgen,ygen)){
 607           fPtVsYGen->Fill(ptgen,ygen);
 608           if(TMath::Abs(ygen)<0.5) fPtVsYGenLimAcc->Fill(ptgen,ygen);
 609           if(isInAcc) fPtVsYGenAcc->Fill(ptgen,ygen);
 610         }
 611         if(TMath::Abs(ygen)<0.9) fPtVsYGenLargeAcc->Fill(ptgen,ygen);
 612       }
 613     }
 614   }
 615 }
 616
 617 //________________________________________________________________________
 618 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){
 619   // Fill histos for candidates with proper charge sign
 620
 621   Bool_t accept=kFALSE;
 622
 623   tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
 624   Double_t pt = tmpRD->Pt();
 625   Double_t minv2 = tmpRD->InvMass2(nProngs,pdgdau);
 626   Double_t mass=TMath::Sqrt(minv2);
 627
 628   if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
 629     Double_t rapid = tmpRD->Y(pdgD);
 630     if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
 631       fMassVsPtVsY->Fill(mass,pt,rapid);
 632       accept=kTRUE;
 633       if(fReadMC){
 634         Int_t signPdg[3]={0,0,0};
 635         for(Int_t iii=0; iii<nProngs; iii++) signPdg[iii]=pdgdau[iii];
 636         Int_t labD = tmpRD->MatchToMC(pdgD,arrayMC,nProngs,signPdg);
 637         if(labD>=0){
 638           AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(TMath::Abs(dgLabels[0])));
 639           if(part){
 640             Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,part,fGoUpToQuark);
 641             if((fPromptFeeddown==kFeeddown && orig==5)|| (fPromptFeeddown==kPrompt && orig==4) || (fPromptFeeddown==kBoth && orig>=4)) {
 642
 643               Int_t pdgCode = TMath::Abs( part->GetPdgCode() );
 644               if(pdgCode==321){
 645                 fMassVsPtVsYSig->Fill(mass,pt,rapid);
 646                 AliAODMCParticle* dmes =  dynamic_cast<AliAODMCParticle*>(arrayMC->At(labD));
 647                 fPtVsYReco->Fill(dmes->Pt(),dmes->Y());
 648               }else{
 649                 fMassVsPtVsYRefl->Fill(mass,pt,rapid);
 650               }
 651             }
 652           }
 653         }else{
 654           fMassVsPtVsYBkg->Fill(mass,pt,rapid);
 655         }
 656       }
 657     }
 658   }
 659
 660   Int_t nRotated=0;
 661   Double_t massRot=0;// calculated later only if candidate is acceptable
 662   Double_t angleProngXY;
 663   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])));
 664   else {
 665     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])));
 666   }
 667   Double_t ptOrig=pt;
 668
 669
 670   Double_t rotStep=(fMaxAngleForRot-fMinAngleForRot)/(fNRotations-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot
 671   Double_t rotStep3=(fMaxAngleForRot3-fMinAngleForRot3)/(fNRotations-1); // -1 is to ensure that the last rotation is done with angle=fMaxAngleForRot
 672
 673   for(Int_t irot=0; irot<fNRotations; irot++){
 674     Double_t phirot=fMinAngleForRot+rotStep*irot;
 675     Double_t tmpx=px[0];
 676     Double_t tmpy=py[0];
 677     Double_t tmpx2=px[2];
 678     Double_t tmpy2=py[2];
 679     px[0]=tmpx*TMath::Cos(phirot)-tmpy*TMath::Sin(phirot);
 680     py[0]=tmpx*TMath::Sin(phirot)+tmpy*TMath::Cos(phirot);
 681     if(pdgD==411 || pdgD==431){
 682       Double_t phirot2=fMinAngleForRot3+rotStep3*irot;
 683       px[2]=tmpx*TMath::Cos(phirot2)-tmpy*TMath::Sin(phirot2);
 684       py[2]=tmpx*TMath::Sin(phirot2)+tmpy*TMath::Cos(phirot2);
 685     }
 686     tmpRD->SetPxPyPzProngs(nProngs,px,py,pz);
 687     pt = tmpRD->Pt();
 688     minv2 = tmpRD->InvMass2(nProngs,pdgdau);
 689     if(minv2>fMinMass*fMinMass && minv2<fMaxMass*fMaxMass){
 690       Double_t rapid = tmpRD->Y(pdgD);
 691       if(fAnalysisCuts->IsInFiducialAcceptance(pt,rapid)){
 692         massRot=TMath::Sqrt(minv2);
 693         fMassVsPtVsYRot->Fill(massRot,pt,rapid);
 694         nRotated++;
 695         fDeltaMass->Fill(massRot-mass);
 696         if(fFullAnalysis){
 697           Double_t pointRot[5]={mass,massRot-mass,ptOrig,pt-ptOrig,angleProngXY};
 698           fDeltaMassFullAnalysis->Fill(pointRot);
 699         }
 700       }
 701     }
 702     px[0]=tmpx;
 703     py[0]=tmpy;
 704     if(pdgD==411 || pdgD==431){
 705       px[2]=tmpx2;
 706       py[2]=tmpy2;
 707     }
 708   }
 709   fNormRotated->Fill(nRotated);
 710
 711   return accept;
 712
 713 }
 714 //________________________________________________________________________
 715 Bool_t AliAnalysisTaskCombinHF::IsTrackSelected(AliAODTrack* track){
 716   // track selection cuts
 717
 718   if(track->Charge()==0) return kFALSE;
 719   if(track->GetID()<0&&!fKeepNegID)return kFALSE;
 720   if(!(track->TestFilterMask(fFilterMask))) return kFALSE;
 721   if(!SelectAODTrack(track,fTrackCutsAll)) return kFALSE;
 722   return kTRUE;
 723 }
 724
 725 //________________________________________________________________________
 726 Bool_t AliAnalysisTaskCombinHF::IsKaon(AliAODTrack* track){
 727   // kaon selection cuts
 728
 729   if(!fPidHF) return kTRUE;
 730   Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);
 731   Double_t mom=track->P();
 732   if(SelectAODTrack(track,fTrackCutsKaon)) {
 733     if(isKaon>=1)    return kTRUE;
 734     if(isKaon<=-1) return kFALSE;
 735     switch(fPIDselCaseZero){// isKaon=0
 736       case 0:
 737       {
 738         return kTRUE;// always accept
 739       }
 740         break;
 741       case 1:
 742       {
 743         if(isKaon>=0 && track->P()>fmaxPforIDKaon)return kTRUE;// accept only if in a compatibility band starting from p=fmaxPforIDKaon
 744       }
 745         break;
 746       case 2:
 747       {
 748         if(track->P()>fmaxPforIDKaon)return kTRUE;
 749         AliPIDResponse *pidResp=fPidHF->GetPidResponse();// more elaborated strategy: asymmetric cuts, with fix momenta and nsigma ranges for the moment
 750         Double_t nsigma=pidResp->NumberOfSigmasTPC(track,AliPID::kKaon);
 751         if(nsigma>-2.&& nsigma<3. && mom<0.6)isKaon=1;
 752         else if(nsigma>-1.&& nsigma<3.&& mom<0.8)isKaon=1;
 753         if(isKaon==1)return kTRUE;
 754       }
 755         break;
 756       default:
 757       {
 758         AliWarning(Form("WRONG CASE OF PID STRATEGY SELECTED: %d (can range from 0 to 2)",fPIDselCaseZero));
 759         return kFALSE;// actually case 0 could be set as the default and return kTRUE
 760       }
 761     }
 762   }
 763
 764   return kFALSE;
 765 }
 766 //_______________________________________________________________________
 767 Bool_t AliAnalysisTaskCombinHF::IsPion(AliAODTrack* track){
 768   // pion selection cuts
 769
 770   if(!fPidHF) return kTRUE;
 771   Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
 772   Double_t mom=track->P();
 773   if(SelectAODTrack(track,fTrackCutsPion)) {
 774     if(isPion>=1)    return kTRUE;
 775     if(isPion<=-1) return kFALSE;
 776     switch(fPIDselCaseZero){// isPion=0
 777       case 0:
 778       {
 779         return kTRUE;// always accept
 780       }
 781         break;
 782       case 1:
 783       {
 784         if(track->P()>fmaxPforIDPion)return kTRUE;// accept only if in a compatibility band starting from p=fmaxPforIDPion
 785       }
 786         break;
 787       case 2:
 788       {
 789         // more elaborated strategy: asymmetric cuts, with fix momenta and nsigma ranges for the moment
 790         if(track->P()>fmaxPforIDPion)return kTRUE;
 791         AliPIDResponse *pidResp=fPidHF->GetPidResponse();
 792         Double_t nsigma=pidResp->NumberOfSigmasTPC(track,AliPID::kPion);
 793         if(nsigma<2.&& nsigma>-3. && mom<0.6)isPion=1;
 794         else if(nsigma<1. && nsigma> -3. && mom<0.8)isPion=1;
 795         if(isPion==1)return kTRUE;
 796       }
 797         break;
 798       default:
 799       {
 800         AliWarning(Form("WRONG CASE OF PID STRATEGY SELECTED: %d (can range from 0 to 2)",fPIDselCaseZero));
 801         return kFALSE;// actually case 0 could be set as the default and return kTRUE
 802       }
 803     }
 804   }
 805
 806   return kFALSE;
 807 }
 808
 809 //________________________________________________________________________
 810 Bool_t AliAnalysisTaskCombinHF::SelectAODTrack(AliAODTrack *track, AliESDtrackCuts *cuts){
 811   // AOD track selection
 812
 813   if(!cuts) return kTRUE;
 814
 815   AliESDtrack esdTrack(track);
 816   // set the TPC cluster info
 817   esdTrack.SetTPCClusterMap(track->GetTPCClusterMap());
 818   esdTrack.SetTPCSharedMap(track->GetTPCSharedMap());
 819   esdTrack.SetTPCPointsF(track->GetTPCNclsF());
 820   if(!cuts->IsSelected(&esdTrack)) return kFALSE;
 821
 822   return kTRUE;
 823 }
 824
 825 //_________________________________________________________________
 826 Bool_t AliAnalysisTaskCombinHF::CheckAcceptance(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
 827   // check if the decay products are in the good eta and pt range
 828   for (Int_t iProng = 0; iProng<nProng; iProng++){
 829     AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
 830     if(!mcPartDaughter) return kFALSE;
 831     Double_t eta = mcPartDaughter->Eta();
 832     Double_t pt = mcPartDaughter->Pt();
 833     if (TMath::Abs(eta) > fEtaAccCut || pt < fPtAccCut) return kFALSE;
 834   }
 835   return kTRUE;
 836 }
 837
 838 //_________________________________________________________________
 839 void AliAnalysisTaskCombinHF::Terminate(Option_t */*option*/)
 840 {
 841   // Terminate analysis
 842   //
 843   if(fDebug > 1) printf("AliAnalysisTaskCombinHF: Terminate() \n");
 844   fOutput = dynamic_cast<TList*> (GetOutputData(1));
 845   if (!fOutput) {
 846     printf("ERROR: fOutput not available\n");
 847     return;
 848   }
 849   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("hNEvents"));
 850   if(fHistNEvents){
 851     printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
 852   }else{
 853     printf("ERROR: fHistNEvents not available\n");
 854     return;
 855   }
 856   return;
 857 }
 858