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