PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, 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 AliAnalysisTaskSEDplus
  20 // AliAnalysisTaskSE for the D+ candidates Invariant Mass Histogram and
  21 //comparison of heavy-flavour decay candidates
  22 // to MC truth (kinematics stored in the AOD)
  23 // Authors: Renu Bala, bala@to.infn.it
  24 // F. Prino, prino@to.infn.it
  25 // G. Ortona, ortona@to.infn.it
  26 /////////////////////////////////////////////////////////////
  27
  28 #include <TClonesArray.h>
  29 #include <TNtuple.h>
  30 #include <TCanvas.h>
  31 #include <TList.h>
  32 #include <TString.h>
  33 #include <TDatabasePDG.h>
  34
  35 #include "AliAnalysisManager.h"
  36 #include "AliRDHFCutsDplustoKpipi.h"
  37 #include "AliAODHandler.h"
  38 #include "AliAODEvent.h"
  39 #include "AliAODVertex.h"
  40 #include "AliAODTrack.h"
  41 #include "AliAODRecoDecayHF3Prong.h"
  42 #include "AliAnalysisVertexingHF.h"
  43 #include "AliAnalysisTaskSE.h"
  44 #include "AliAnalysisTaskSEDplus.h"
  45 #include "AliNormalizationCounter.h"
  46
  47 ClassImp(AliAnalysisTaskSEDplus)
  48
  49
  50 //________________________________________________________________________
  51 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
  52 AliAnalysisTaskSE(),
  53   fOutput(0),
  54   fHistNEvents(0),
  55   fPtVsMass(0),
  56   fPtVsMassTC(0),
  57   fYVsPt(0),
  58   fYVsPtTC(0),
  59   fYVsPtSig(0),
  60   fYVsPtSigTC(0),
  61   fNtupleDplus(0),
  62   fUpmasslimit(1.965),
  63   fLowmasslimit(1.765),
  64   fNPtBins(0),
  65   fBinWidth(0.002),
  66   fListCuts(0),
  67   fRDCutsProduction(0),
  68   fRDCutsAnalysis(0),
  69   fCounter(0),
  70   fFillNtuple(kFALSE),
  71   fReadMC(kFALSE),
  72   fUseStrangeness(kFALSE),
  73   fUseBit(kTRUE),
  74   fCutsDistr(kFALSE),
  75   fDoImpPar(kFALSE),
  76   fNImpParBins(400),
  77   fLowerImpPar(-2000.),
  78   fHigherImpPar(2000.),
  79   fDoLS(0)
  80 {
  81    // Default constructor
  82 }
  83
  84 //________________________________________________________________________
  85 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,AliRDHFCutsDplustoKpipi *dpluscutsana,AliRDHFCutsDplustoKpipi *dpluscutsprod,Bool_t fillNtuple):
  86   AliAnalysisTaskSE(name),
  87   fOutput(0),
  88   fHistNEvents(0),
  89   fPtVsMass(0),
  90   fPtVsMassTC(0),
  91   fYVsPt(0),
  92   fYVsPtTC(0),
  93   fYVsPtSig(0),
  94   fYVsPtSigTC(0),
  95   fNtupleDplus(0),
  96   fUpmasslimit(1.965),
  97   fLowmasslimit(1.765),
  98   fNPtBins(0),
  99   fBinWidth(0.002),
 100   fListCuts(0),
 101   fRDCutsProduction(dpluscutsprod),
 102   fRDCutsAnalysis(dpluscutsana),
 103   fCounter(0),
 104   fFillNtuple(fillNtuple),
 105   fReadMC(kFALSE),
 106   fUseStrangeness(kFALSE),
 107   fUseBit(kTRUE),
 108   fCutsDistr(kFALSE),
 109   fDoImpPar(kFALSE),
 110   fNImpParBins(400),
 111   fLowerImpPar(-2000.),
 112   fHigherImpPar(2000.),
 113   fDoLS(0)
 114 {
 115   //
 116   // Standrd constructor
 117   //
 118   fNPtBins=fRDCutsAnalysis->GetNPtBins();
 119   // Default constructor
 120    // Output slot #1 writes into a TList container
 121   DefineOutput(1,TList::Class());  //My private output
 122  // Output slot #2 writes cut to private output
 123   //  DefineOutput(2,AliRDHFCutsDplustoKpipi::Class());
 124   DefineOutput(2,TList::Class());
 125 // Output slot #3 writes cut to private output
 126   DefineOutput(3,AliNormalizationCounter::Class());
 127
 128   if(fFillNtuple){
 129     // Output slot #4 writes into a TNtuple container
 130     DefineOutput(4,TNtuple::Class());  //My private output
 131   }
 132 }
 133
 134 //________________________________________________________________________
 135 AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
 136 {
 137   //
 138   // Destructor
 139   //
 140   if (fOutput) {
 141     delete fOutput;
 142     fOutput = 0;
 143   }
 144   if(fHistNEvents){
 145     delete fHistNEvents;
 146     fHistNEvents=0;
 147   }
 148   for(Int_t i=0;i<3;i++){
 149     if(fHistCentrality[i]){delete fHistCentrality[i]; fHistCentrality[i]=0;}
 150   }
 151
 152   for(Int_t i=0;i<3*fNPtBins;i++){
 153     if(fMassHist[i]){ delete fMassHist[i]; fMassHist[i]=0;}
 154     if(fCosPHist[i]){ delete fCosPHist[i]; fCosPHist[i]=0;}
 155     if(fDLenHist[i]){ delete fDLenHist[i]; fDLenHist[i]=0;}
 156     if(fSumd02Hist[i]){ delete fSumd02Hist[i]; fSumd02Hist[i]=0;}
 157     if(fSigVertHist[i]){ delete fSigVertHist[i]; fSigVertHist[i]=0;}
 158     if(fPtMaxHist[i]){ delete fPtMaxHist[i]; fPtMaxHist[i]=0;}
 159     if(fPtKHist[i]){ delete fPtKHist[i]; fPtKHist[i]=0;}
 160     if(fPtpi1Hist[i]){ delete fPtpi1Hist[i]; fPtpi1Hist[i]=0;}
 161     if(fPtpi2Hist[i]){ delete fPtpi2Hist[i]; fPtpi2Hist[i]=0;}
 162     if(fDCAHist[i]){ delete fDCAHist[i]; fDCAHist[i]=0;}
 163     if(fMassHistTC[i]){ delete fMassHistTC[i]; fMassHistTC[i]=0;}
 164     if(fMassHistTCPlus[i]){ delete fMassHistTCPlus[i]; fMassHistTCPlus[i]=0;}
 165     if(fMassHistTCMinus[i]){ delete fMassHistTCMinus[i]; fMassHistTCMinus[i]=0;}
 166
 167     if(fDLxy[i]){delete fDLxy[i]; fDLxy[i]=0;}
 168     if(fDLxyTC[i]){delete fDLxyTC[i]; fDLxyTC[i]=0;}
 169     if(fCosxy[i]){delete fCosxy[i]; fCosxy[i]=0;}
 170     if(fCosxyTC[i]){delete fCosxyTC[i]; fCosxyTC[i]=0;}
 171     if(fMassHistLS[i]){ delete fMassHistLS[i]; fMassHistLS[i]=0;}
 172     if(fCosPHistLS[i]){ delete fCosPHistLS[i]; fCosPHistLS[i]=0;}
 173     if(fDLenHistLS[i]){ delete fDLenHistLS[i]; fDLenHistLS[i]=0;}
 174     if(fSumd02HistLS[i]){ delete fSumd02HistLS[i]; fSumd02HistLS[i]=0;}
 175     if(fSigVertHistLS[i]){ delete fSigVertHistLS[i]; fSigVertHistLS[i]=0;}
 176     if(fPtMaxHistLS[i]){ delete fPtMaxHistLS[i]; fPtMaxHistLS[i]=0;}
 177     if(fDCAHistLS[i]){ delete fDCAHistLS[i]; fDCAHistLS[i]=0;}
 178     if(fMassHistLSTC[i]){ delete fMassHistLSTC[i]; fMassHistLSTC[i]=0;}
 179   }
 180
 181   for(Int_t i=0;i<3;i++){
 182     if(fCorreld0Kd0pi[i]){ delete fCorreld0Kd0pi[i]; fCorreld0Kd0pi[i]=0;}
 183   }
 184
 185   if(fPtVsMass){
 186     delete fPtVsMass;
 187     fPtVsMass=0;
 188   }
 189   if(fPtVsMassTC){
 190     delete fPtVsMassTC;
 191     fPtVsMassTC=0;
 192   }
 193   if(fYVsPt){
 194     delete fYVsPt;
 195     fYVsPt=0;
 196   }
 197   if(fYVsPtTC){
 198     delete fYVsPtTC;
 199     fYVsPtTC=0;
 200   }
 201   if(fYVsPtSig){
 202     delete fYVsPtSig;
 203     fYVsPtSig=0;
 204   }
 205   if(fYVsPtSigTC){
 206     delete fYVsPtSigTC;
 207     fYVsPtSigTC=0;
 208   }
 209   if(fNtupleDplus){
 210     delete fNtupleDplus;
 211     fNtupleDplus=0;
 212   }
 213   if (fListCuts) {
 214     delete fListCuts;
 215     fListCuts = 0;
 216   }
 217   if(fRDCutsProduction){
 218     delete fRDCutsProduction;
 219     fRDCutsProduction = 0;
 220   }
 221   if(fRDCutsAnalysis){
 222     delete fRDCutsAnalysis;
 223     fRDCutsAnalysis = 0;
 224   }
 225   for(Int_t i=0; i<5; i++){
 226     delete fHistMassPtImpParTC[i];
 227   }
 228   if(fCounter){
 229     delete fCounter;
 230     fCounter = 0;
 231   }
 232
 233
 234 }
 235 //_________________________________________________________________
 236 void  AliAnalysisTaskSEDplus::SetMassLimits(Float_t range){
 237   // set invariant mass limits
 238   Float_t bw=GetBinWidth();
 239   fUpmasslimit = 1.865+range;
 240   fLowmasslimit = 1.865-range;
 241   SetBinWidth(bw);
 242 }
 243 //_________________________________________________________________
 244 void  AliAnalysisTaskSEDplus::SetMassLimits(Float_t lowlimit, Float_t uplimit){
 245   // set invariant mass limits
 246   if(uplimit>lowlimit)
 247     {
 248       Float_t bw=GetBinWidth();
 249       fUpmasslimit = lowlimit;
 250       fLowmasslimit = uplimit;
 251       SetBinWidth(bw);
 252     }
 253 }
 254 //________________________________________________________________
 255 void AliAnalysisTaskSEDplus::SetBinWidth(Float_t w){
 256   Float_t width=w;
 257   Int_t nbins=(Int_t)((fUpmasslimit-fLowmasslimit)/width+0.5);
 258   Int_t missingbins=4-nbins%4;
 259   nbins=nbins+missingbins;
 260   width=(fUpmasslimit-fLowmasslimit)/nbins;
 261   if(missingbins!=0){
 262     printf("AliAnalysisTaskSEDplus::SetBinWidth: W-bin width of %f will produce histograms not rebinnable by 4. New width set to %f\n",w,width);
 263   }
 264   else{
 265     if(fDebug>1) printf("AliAnalysisTaskSEDplus::SetBinWidth: width set to %f\n",width);
 266   }
 267   fBinWidth=width;
 268 }
 269 //_________________________________________________________________
 270 Int_t AliAnalysisTaskSEDplus::GetNBinsHistos(){
 271   return (Int_t)((fUpmasslimit-fLowmasslimit)/fBinWidth+0.5);
 272 }
 273 //_________________________________________________________________
 274 void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS){
 275   //
 276   //
 277   // Fill the Like Sign histograms
 278   //
 279   if(fDebug>1)printf("started LS\n");
 280
 281   //histograms for like sign
 282   Int_t nbins=GetNBinsHistos();;
 283   TH1F *histLSPlus = new TH1F("LSPlus","LSPlus",nbins,fLowmasslimit,fUpmasslimit);
 284   TH1F *histLSMinus = new TH1F("LSMinus","LSMinus",nbins,fLowmasslimit,fUpmasslimit);
 285
 286   Int_t nPosTrks=0,nNegTrks=0;
 287   Int_t nOStriplets = arrayOppositeSign->GetEntriesFast();
 288   Int_t nDplusLS=0;
 289   Int_t nDminusLS=0;
 290   Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
 291   Int_t index=0;
 292
 293   // loop over like sign candidates
 294   for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
 295     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
 296     if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts))continue;
 297     Bool_t unsetvtx=kFALSE;
 298     if(!d->GetOwnPrimaryVtx()) {
 299       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 300       unsetvtx=kTRUE;
 301     }
 302     Double_t ptCand = d->Pt();
 303     Int_t iPtBin = fRDCutsAnalysis->PtBin(ptCand);
 304     if(iPtBin<0)continue;
 305     Int_t sign= d->GetCharge();
 306     if(sign>0){
 307       nPosTrks++;
 308     }else{
 309       nNegTrks++;
 310     }
 311     //    if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)){
 312     fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
 313     Int_t passTightCuts=fRDCutsAnalysis->GetIsSelectedCuts();
 314     if(passTightCuts>0){
 315       Float_t invMass = d->InvMassDplus();
 316       index=GetLSHistoIndex(iPtBin);
 317       fMassHistLS[index+1]->Fill(invMass);
 318       if(sign>0){
 319         histLSPlus->Fill(invMass);
 320         nDplusLS++;
 321       }else{
 322         histLSMinus->Fill(invMass);
 323         nDminusLS++;
 324       }
 325       if(fCutsDistr){
 326         Double_t dlen=d->DecayLength();
 327         Double_t cosp=d->CosPointingAngle();
 328         Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 329         Double_t dca=d->GetDCA();
 330         Double_t sigvert=d->GetSigmaVert();
 331         Double_t ptmax=0;
 332         for(Int_t i=0;i<3;i++){
 333           if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
 334         }
 335         fCosPHistLS[iPtBin]->Fill(cosp);
 336         fDLenHistLS[iPtBin]->Fill(dlen);
 337         fSumd02HistLS[iPtBin]->Fill(sumD02);
 338         fSigVertHistLS[iPtBin]->Fill(sigvert);
 339         fPtMaxHistLS[iPtBin]->Fill(ptmax);
 340         fDCAHistLS[iPtBin]->Fill(dca);
 341       }
 342     }
 343     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 344   }
 345   //wei:normalized to the number of combinations (production)
 346   //wei2:normalized to the number of  LS/OS (production)
 347   //wei3:normalized to the number of  LS/OS (analysis)
 348   //wei4:normalized to the number of combinations (analysis)
 349   Float_t wei2=0;
 350   if(nLikeSign!=0)wei2 = (Float_t)nOStriplets/(Float_t)nLikeSign;
 351   Float_t wei3=0;
 352   if(nDplusLS!=0)wei3 = (Float_t)nDplusOS/(Float_t)(nDplusLS+nDminusLS);
 353   Float_t weiplus=1.,weiminus=1.;
 354   Float_t wei4plus=1.,wei4minus=1.;
 355   //wei* should be automatically protected, since to get a triplet there must be at least 3 good tracks in the event
 356   if(nPosTrks>2)weiplus=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.);
 357   if(nDplusLS>2)wei4plus=3.*(Float_t)nDminusLS/((Float_t)nDplusLS-2.);
 358   if(nNegTrks>2)weiminus=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.);
 359   if(nDminusLS>2)wei4minus=3.*(Float_t)nDplusLS/((Float_t)nDminusLS-2.);
 360
 361   fMassHistLS[index]->Add(histLSPlus,weiplus);
 362   fMassHistLS[index]->Add(histLSMinus,weiminus);
 363   fMassHistLS[index+2]->Add(histLSPlus,wei2);
 364   fMassHistLS[index+2]->Add(histLSMinus,wei2);
 365   fMassHistLS[index+3]->Add(histLSPlus,wei3);
 366   fMassHistLS[index+3]->Add(histLSMinus,wei3);
 367   fMassHistLS[index+4]->Add(histLSPlus,wei4plus);
 368   fMassHistLS[index+4]->Add(histLSMinus,wei4minus);
 369
 370   delete histLSPlus;histLSPlus=0;
 371   delete histLSMinus;histLSMinus=0;
 372
 373   if(fDebug>1) printf("LS analysis terminated\n");
 374 }
 375
 376 //__________________________________________
 377 void AliAnalysisTaskSEDplus::Init(){
 378   //
 379   // Initialization
 380   //
 381   if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
 382
 383   //PostData(2,fRDCutsloose);//we should then put those cuts in a tlist if we have more than 1
 384   fListCuts=new TList();
 385   AliRDHFCutsDplustoKpipi *production = new AliRDHFCutsDplustoKpipi(*fRDCutsProduction);
 386   production->SetName("ProductionCuts");
 387   AliRDHFCutsDplustoKpipi *analysis = new AliRDHFCutsDplustoKpipi(*fRDCutsAnalysis);
 388   analysis->SetName("AnalysisCuts");
 389
 390   fListCuts->Add(production);
 391   fListCuts->Add(analysis);
 392   PostData(2,fListCuts);
 393
 394   return;
 395 }
 396
 397 //________________________________________________________________________
 398 void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
 399 {
 400   // Create the output container
 401   //
 402   if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");
 403
 404   // Several histograms are more conveniently managed in a TList
 405   fOutput = new TList();
 406   fOutput->SetOwner();
 407   fOutput->SetName("OutputHistos");
 408
 409   TString hisname;
 410   Int_t index=0;
 411   Int_t nbins=GetNBinsHistos();
 412   fHistCentrality[0]=new TH1F("centrality","centrality",100,0.5,30000.5);
 413   fHistCentrality[1]=new TH1F("centrality(selectedCent)","centrality(selectedCent)",100,0.5,30000.5);
 414   fHistCentrality[2]=new TH1F("centrality(OutofCent)","centrality(OutofCent)",100,0.5,30000.5);
 415   for(Int_t i=0;i<3;i++){
 416     fHistCentrality[i]->Sumw2();
 417     fOutput->Add(fHistCentrality[i]);
 418   }
 419   for(Int_t i=0;i<fNPtBins;i++){
 420
 421     index=GetHistoIndex(i);
 422
 423     hisname.Form("hMassPt%d",i);
 424     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 425     fMassHist[index]->Sumw2();
 426     hisname.Form("hCosPAllPt%d",i);
 427     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
 428     fCosPHist[index]->Sumw2();
 429     hisname.Form("hDLenAllPt%d",i);
 430     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
 431     fDLenHist[index]->Sumw2();
 432     hisname.Form("hSumd02AllPt%d",i);
 433     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
 434     fSumd02Hist[index]->Sumw2();
 435     hisname.Form("hSigVertAllPt%d",i);
 436     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 437     fSigVertHist[index]->Sumw2();
 438     hisname.Form("hPtMaxAllPt%d",i);
 439     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 440     fPtMaxHist[index]->Sumw2();
 441     hisname.Form("hPtKPt%d",i);
 442     fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 443     fPtKHist[index]->Sumw2();
 444     hisname.Form("hPtpi1Pt%d",i);
 445     fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 446     fPtpi1Hist[index]->Sumw2();
 447     hisname.Form("hPtpi2Pt%d",i);
 448     fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 449     fPtpi2Hist[index]->Sumw2();
 450     hisname.Form("hDCAAllPt%d",i);
 451     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 452     fDCAHist[index]->Sumw2();
 453
 454     hisname.Form("hDLxyPt%d",i);
 455     fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
 456     fDLxy[index]->Sumw2();
 457     hisname.Form("hCosxyPt%d",i);
 458     fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
 459     fCosxy[index]->Sumw2();
 460     hisname.Form("hDLxyPt%dTC",i);
 461     fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
 462     fDLxyTC[index]->Sumw2();
 463     hisname.Form("hCosxyPt%dTC",i);
 464     fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
 465     fCosxyTC[index]->Sumw2();
 466
 467     hisname.Form("hMassPt%dTC",i);
 468     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 469     fMassHistTC[index]->Sumw2();
 470     hisname.Form("hMassPt%dTCPlus",i);
 471     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 472     fMassHistTCPlus[index]->Sumw2();
 473     hisname.Form("hMassPt%dTCMinus",i);
 474     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 475     fMassHistTCMinus[index]->Sumw2();
 476
 477
 478
 479     index=GetSignalHistoIndex(i);
 480     hisname.Form("hSigPt%d",i);
 481     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 482     fMassHist[index]->Sumw2();
 483     hisname.Form("hCosPSigPt%d",i);
 484     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
 485     fCosPHist[index]->Sumw2();
 486     hisname.Form("hDLenSigPt%d",i);
 487     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
 488     fDLenHist[index]->Sumw2();
 489     hisname.Form("hSumd02SigPt%d",i);
 490     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
 491     fSumd02Hist[index]->Sumw2();
 492     hisname.Form("hSigVertSigPt%d",i);
 493     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 494     fSigVertHist[index]->Sumw2();
 495     hisname.Form("hPtMaxSigPt%d",i);
 496     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 497     fPtMaxHist[index]->Sumw2();
 498     hisname.Form("hPtKSigPt%d",i);
 499     fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 500     fPtKHist[index]->Sumw2();
 501     hisname.Form("hPtpi1SigPt%d",i);
 502     fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 503     fPtpi1Hist[index]->Sumw2();
 504     hisname.Form("hPtpi2SigPt%d",i);
 505     fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 506     fPtpi2Hist[index]->Sumw2();
 507
 508     hisname.Form("hDCASigPt%d",i);
 509     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 510     fDCAHist[index]->Sumw2();
 511
 512     hisname.Form("hDLxySigPt%d",i);
 513     fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
 514     fDLxy[index]->Sumw2();
 515     hisname.Form("hCosxySigPt%d",i);
 516     fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
 517     fCosxy[index]->Sumw2();
 518     hisname.Form("hDLxySigPt%dTC",i);
 519     fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
 520     fDLxyTC[index]->Sumw2();
 521     hisname.Form("hCosxySigPt%dTC",i);
 522     fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
 523     fCosxyTC[index]->Sumw2();
 524     hisname.Form("hSigPt%dTC",i);
 525     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 526     fMassHistTC[index]->Sumw2();
 527     hisname.Form("hSigPt%dTCPlus",i);
 528     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 529     fMassHistTCPlus[index]->Sumw2();
 530     hisname.Form("hSigPt%dTCMinus",i);
 531     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 532     fMassHistTCMinus[index]->Sumw2();
 533
 534
 535     index=GetBackgroundHistoIndex(i);
 536     hisname.Form("hBkgPt%d",i);
 537     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 538     fMassHist[index]->Sumw2();
 539     hisname.Form("hCosPBkgPt%d",i);
 540     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
 541     fCosPHist[index]->Sumw2();
 542     hisname.Form("hDLenBkgPt%d",i);
 543     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
 544     fDLenHist[index]->Sumw2();
 545     hisname.Form("hSumd02BkgPt%d",i);
 546     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
 547     fSumd02Hist[index]->Sumw2();
 548     hisname.Form("hSigVertBkgPt%d",i);
 549     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 550     fSigVertHist[index]->Sumw2();
 551     hisname.Form("hPtMaxBkgPt%d",i);
 552     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 553     fPtMaxHist[index]->Sumw2();
 554     hisname.Form("hPtKBkgPt%d",i);
 555     fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 556     fPtKHist[index]->Sumw2();
 557     hisname.Form("hPtpi1BkgPt%d",i);
 558     fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 559     fPtpi1Hist[index]->Sumw2();
 560     hisname.Form("hPtpi2BkgPt%d",i);
 561     fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 562     fPtpi2Hist[index]->Sumw2();
 563     hisname.Form("hDCABkgPt%d",i);
 564     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 565     fDCAHist[index]->Sumw2();
 566
 567     hisname.Form("hDLxyBkgPt%d",i);
 568     fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
 569     fDLxy[index]->Sumw2();
 570     hisname.Form("hCosxyBkgPt%d",i);
 571     fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
 572     fCosxy[index]->Sumw2();
 573     hisname.Form("hDLxyBkgPt%dTC",i);
 574     fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,10.);
 575     fDLxyTC[index]->Sumw2();
 576     hisname.Form("hCosxyBkgPt%dTC",i);
 577     fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.85,1.);
 578     fCosxyTC[index]->Sumw2();
 579
 580
 581     hisname.Form("hBkgPt%dTC",i);
 582     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 583     fMassHistTC[index]->Sumw2();
 584     hisname.Form("hBkgPt%dTCPlus",i);
 585     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 586     fMassHistTCPlus[index]->Sumw2();
 587     hisname.Form("hBkgPt%dTCMinus",i);
 588     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 589     fMassHistTCMinus[index]->Sumw2();
 590   }
 591
 592
 593   for(Int_t i=0; i<3*fNPtBins; i++){
 594     fOutput->Add(fMassHist[i]);
 595     if(fCutsDistr){
 596       fOutput->Add(fCosPHist[i]);
 597       fOutput->Add(fDLenHist[i]);
 598       fOutput->Add(fSumd02Hist[i]);
 599       fOutput->Add(fSigVertHist[i]);
 600       fOutput->Add(fPtMaxHist[i]);
 601       fOutput->Add(fPtKHist[i]);
 602       fOutput->Add(fPtpi1Hist[i]);
 603       fOutput->Add(fPtpi2Hist[i]);
 604       fOutput->Add(fDCAHist[i]);
 605       fOutput->Add(fDLxy[i]);
 606       fOutput->Add(fDLxyTC[i]);
 607       fOutput->Add(fCosxy[i]);
 608       fOutput->Add(fCosxyTC[i]);
 609     }
 610     fOutput->Add(fMassHistTC[i]);
 611     fOutput->Add(fMassHistTCPlus[i]);
 612     fOutput->Add(fMassHistTCMinus[i]);
 613
 614   }
 615
 616   if(fCutsDistr){
 617     fCorreld0Kd0pi[0]=new TH2F("hCorreld0Kd0piAll","",100,-0.02,0.02,100,-0.02,0.02);
 618     fCorreld0Kd0pi[1]=new TH2F("hCorreld0Kd0piSig","",100,-0.02,0.02,100,-0.02,0.02);
 619     fCorreld0Kd0pi[2]=new TH2F("hCorreld0Kd0piBkg","",100,-0.02,0.02,100,-0.02,0.02);
 620     for(Int_t i=0; i<3; i++){
 621       fCorreld0Kd0pi[i]->Sumw2();
 622       fOutput->Add(fCorreld0Kd0pi[i]);
 623     }
 624   }
 625
 626   fHistNEvents = new TH1F("fHistNEvents", "number of events ",9,-0.5,8.5);
 627   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 628   fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with good vertex");
 629   fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents with PbPb HM trigger");
 630   fHistNEvents->GetXaxis()->SetBinLabel(4,"no. of Rejected pileup events");
 631   fHistNEvents->GetXaxis()->SetBinLabel(5,"no. of candidate");
 632   fHistNEvents->GetXaxis()->SetBinLabel(6,"no. of D+ after loose cuts");
 633   fHistNEvents->GetXaxis()->SetBinLabel(7,"no. of D+ after tight cuts");
 634   fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events");
 635   fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
 636   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
 637   fHistNEvents->Sumw2();
 638   fHistNEvents->SetMinimum(0);
 639   fOutput->Add(fHistNEvents);
 640
 641   fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
 642   fPtVsMassTC=new TH2F("hPtVsMassTC","PtVsMass (analysis cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
 643   fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
 644   fYVsPtTC=new TH2F("hYVsPtTC","YvsPt (analysis cuts)",40,0.,20.,80,-2.,2.);
 645   fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
 646   fYVsPtSigTC=new TH2F("hYVsPtSigTC","YvsPt (MC, only Sig, analysis cuts)",40,0.,20.,80,-2.,2.);
 647
 648   fOutput->Add(fPtVsMass);
 649   fOutput->Add(fPtVsMassTC);
 650   fOutput->Add(fYVsPt);
 651   fOutput->Add(fYVsPtTC);
 652   fOutput->Add(fYVsPtSig);
 653   fOutput->Add(fYVsPtSigTC);
 654
 655
 656   //Counter for Normalization
 657   TString normName="NormalizationCounter";
 658   AliAnalysisDataContainer *cont = GetOutputSlot(3)->GetContainer();
 659   if(cont)normName=(TString)cont->GetName();
 660   fCounter = new AliNormalizationCounter(normName.Data());
 661   fCounter->Init();
 662
 663   if(fDoLS) CreateLikeSignHistos();
 664   if(fDoImpPar) CreateImpactParameterHistos();
 665
 666   if(fFillNtuple){
 667     OpenFile(4); // 4 is the slot number of the ntuple
 668
 669     fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:PtTrue:VxTrue:VyTrue:VzTrue:Ptpi:PtK:Ptpi2:PtRec:PointingAngle:DecLeng:VxRec:VyRec:VzRec:InvMass:sigvert:d0Pi:d0K:d0Pi2:dca:d0square");
 670
 671   }
 672
 673   return;
 674 }
 675
 676 //________________________________________________________________________
 677 void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
 678 {
 679   // Execute analysis for current event:
 680   // heavy flavor candidates association to MC truth
 681
 682   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 683
 684
 685
 686   TClonesArray *array3Prong = 0;
 687   TClonesArray *arrayLikeSign =0;
 688   if(!aod && AODEvent() && IsStandardAOD()) {
 689     // In case there is an AOD handler writing a standard AOD, use the AOD
 690     // event in memory rather than the input (ESD) event.
 691     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 692      // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 693      // have to taken from the AOD event hold by the AliAODExtension
 694     AliAODHandler* aodHandler = (AliAODHandler*)
 695       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 696     if(aodHandler->GetExtensions()) {
 697       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 698       AliAODEvent *aodFromExt = ext->GetAOD();
 699       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 700       arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
 701     }
 702   } else if(aod) {
 703     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 704     arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
 705   }
 706
 707   if(!aod || !array3Prong) {
 708     printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
 709     return;
 710   }
 711   if(!arrayLikeSign && fDoLS) {
 712     printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
 713     return;
 714   }
 715
 716
 717   // fix for temporary bug in ESDfilter
 718   // the AODs with null vertex pointer didn't pass the PhysSel
 719   if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
 720   fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC);
 721   fHistNEvents->Fill(0); // count event
 722
 723   Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
 724   Float_t centrality=aod->GetNTracks();//fRDCutsAnalysis->GetCentrality(aod);
 725   fHistCentrality[0]->Fill(centrality);
 726   // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
 727   TString trigclass=aod->GetFiredTriggerClasses();
 728   if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")||trigclass.Contains("C0SMH-B-NOPF-ALL")) fHistNEvents->Fill(2);
 729   if(fRDCutsAnalysis->GetWhyRejection()==1)fHistNEvents->Fill(3);
 730   if(fRDCutsAnalysis->GetWhyRejection()==2){fHistNEvents->Fill(7);fHistCentrality[2]->Fill(centrality);}
 731
 732   // Post the data already here
 733   PostData(1,fOutput);
 734   if(!isEvSel)return;
 735
 736   fHistCentrality[1]->Fill(centrality);
 737   fHistNEvents->Fill(1);
 738
 739   TClonesArray *arrayMC=0;
 740   AliAODMCHeader *mcHeader=0;
 741
 742   // AOD primary vertex
 743   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 744   //    vtx1->Print();
 745    TString primTitle = vtx1->GetTitle();
 746    //if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0)fHistNEvents->Fill(2);
 747
 748   // load MC particles
 749   if(fReadMC){
 750
 751     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 752     if(!arrayMC) {
 753       printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
 754       return;
 755     }
 756
 757   // load MC header
 758     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 759     if(!mcHeader) {
 760     printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
 761     return;
 762     }
 763   }
 764
 765   Int_t n3Prong = array3Prong->GetEntriesFast();
 766   //  printf("Number of D+->Kpipi: %d and of tracks: %d\n",n3Prong,aod->GetNumberOfTracks());
 767
 768
 769   Int_t nOS=0;
 770   Int_t index;
 771   Int_t pdgDgDplustoKpipi[3]={321,211,211};
 772
 773   if(fDoLS>1){
 774     for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 775       AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 776       if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
 777         if(fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod))nOS++;
 778       }
 779     }
 780   }else{
 781   // Double_t cutsDplus[12]={0.2,0.4,0.4,0.,0.,0.01,0.06,0.02,0.,0.85,0.,10000000000.};//TO REMOVE
 782   //Double_t *cutsDplus = new (Double_t*)fRDCuts->GetCuts();
 783   Int_t nSelectedloose=0,nSelectedtight=0;
 784   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 785     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 786     fHistNEvents->Fill(4);
 787     if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
 788       fHistNEvents->Fill(8);
 789       continue;
 790     }
 791     Bool_t unsetvtx=kFALSE;
 792     if(!d->GetOwnPrimaryVtx()){
 793       d->SetOwnPrimaryVtx(vtx1);
 794       unsetvtx=kTRUE;
 795     }
 796
 797     if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)) {
 798
 799       Double_t ptCand = d->Pt();
 800       Int_t iPtBin = fRDCutsProduction->PtBin(ptCand);
 801
 802       Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
 803       Bool_t recVtx=kFALSE;
 804       AliAODVertex *origownvtx=0x0;
 805       if(fRDCutsProduction->GetIsPrimaryWithoutDaughters()){
 806         if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
 807         if(fRDCutsProduction->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
 808         else fRDCutsProduction->CleanOwnPrimaryVtx(d,aod,origownvtx);
 809       }
 810
 811       Int_t labDp=-1;
 812       Bool_t isPrimary=kTRUE;
 813       Float_t deltaPx=0.;
 814       Float_t deltaPy=0.;
 815       Float_t deltaPz=0.;
 816       Float_t truePt=0.;
 817       Float_t xDecay=0.;
 818       Float_t yDecay=0.;
 819       Float_t zDecay=0.;
 820       Float_t pdgCode=-2;
 821       Float_t trueImpParXY=0.;
 822       if(fReadMC){
 823         labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
 824         if(labDp>=0){
 825           AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
 826           if(CheckOrigin(arrayMC,partDp)==5) isPrimary=kFALSE;
 827           AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
 828           deltaPx=partDp->Px()-d->Px();
 829           deltaPy=partDp->Py()-d->Py();
 830           deltaPz=partDp->Pz()-d->Pz();
 831           truePt=partDp->Pt();
 832           xDecay=dg0->Xv();
 833           yDecay=dg0->Yv();
 834           zDecay=dg0->Zv();
 835           pdgCode=TMath::Abs(partDp->GetPdgCode());
 836           if(!isPrimary){
 837             trueImpParXY=GetTrueImpactParameter(mcHeader,arrayMC,partDp)*10000.;
 838           }
 839         }else{
 840           pdgCode=-1;
 841         }
 842       }
 843
 844       Double_t invMass=d->InvMassDplus();
 845       Double_t rapid=d->YDplus();
 846       fYVsPt->Fill(ptCand,rapid);
 847       if(passTightCuts) {fYVsPtTC->Fill(ptCand,rapid);nOS++;}
 848       Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
 849       if(isFidAcc){
 850         fPtVsMass->Fill(invMass,ptCand);
 851         if(passTightCuts) fPtVsMassTC->Fill(invMass,ptCand);
 852       }
 853       Float_t tmp[24];
 854       Double_t  dlen=d->DecayLength();
 855       Double_t cosp=d->CosPointingAngle();
 856       Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 857       Double_t dca=d->GetDCA();
 858       Double_t sigvert=d->GetSigmaVert();
 859       Double_t ptmax=0;
 860       for(Int_t i=0;i<3;i++){
 861         if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
 862       }
 863       Double_t impparXY=d->ImpParXY()*10000.;
 864       Double_t arrayForSparse[3]={invMass,ptCand,impparXY};
 865       Double_t arrayForSparseTrue[3]={invMass,ptCand,trueImpParXY};
 866       if(fFillNtuple){
 867         tmp[0]=pdgCode;
 868         tmp[1]=deltaPx;
 869         tmp[2]=deltaPy;
 870         tmp[3]=deltaPz;
 871         tmp[4]=truePt;
 872         tmp[5]=xDecay;
 873         tmp[6]=yDecay;
 874         tmp[7]=zDecay;
 875         tmp[8]=d->PtProng(0);
 876         tmp[9]=d->PtProng(1);
 877         tmp[10]=d->PtProng(2);
 878         tmp[11]=d->Pt();
 879         tmp[12]=cosp;
 880         tmp[13]=dlen;
 881         tmp[14]=d->Xv();
 882         tmp[15]=d->Yv();
 883         tmp[16]=d->Zv();
 884         tmp[17]=d->InvMassDplus();
 885         tmp[18]=sigvert;
 886         tmp[19]=d->Getd0Prong(0);
 887         tmp[20]=d->Getd0Prong(1);
 888         tmp[21]=d->Getd0Prong(2);
 889         tmp[22]=dca;
 890         tmp[23]=d->Prodd0d0();
 891         fNtupleDplus->Fill(tmp);
 892         PostData(4,fNtupleDplus);
 893       }
 894       if(iPtBin>=0){
 895         Float_t dlxy=d->NormalizedDecayLengthXY();
 896         Float_t cxy=d->CosPointingAngleXY();
 897         index=GetHistoIndex(iPtBin);
 898         if(isFidAcc){
 899           fHistNEvents->Fill(5);
 900           nSelectedloose++;
 901           fMassHist[index]->Fill(invMass);
 902           if(fCutsDistr){
 903             fCosPHist[index]->Fill(cosp);
 904             fDLenHist[index]->Fill(dlen);
 905             fSumd02Hist[index]->Fill(sumD02);
 906             fSigVertHist[index]->Fill(sigvert);
 907             fPtMaxHist[index]->Fill(ptmax);
 908             fPtKHist[index]->Fill(d->PtProng(1));
 909             fPtpi1Hist[index]->Fill(d->PtProng(0));
 910             fPtpi2Hist[index]->Fill(d->PtProng(2));
 911             fDCAHist[index]->Fill(dca);
 912             fDLxy[index]->Fill(dlxy);
 913             fCosxy[index]->Fill(cxy);
 914             fCorreld0Kd0pi[0]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
 915                                     d->Getd0Prong(2)*d->Getd0Prong(1));
 916           }
 917           if(passTightCuts){ fHistNEvents->Fill(6);
 918             nSelectedtight++;
 919             fMassHistTC[index]->Fill(invMass);
 920             if(fCutsDistr){
 921               fDLxyTC[index]->Fill(dlxy);
 922               fCosxyTC[index]->Fill(cxy);
 923             }
 924             if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
 925             else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
 926             if(fDoImpPar){
 927               fHistMassPtImpParTC[0]->Fill(arrayForSparse);
 928             }
 929           }
 930         }
 931
 932         if(fReadMC){
 933           //  if(fCutsDistr){
 934           if(labDp>=0) {
 935             index=GetSignalHistoIndex(iPtBin);
 936             if(isFidAcc){
 937               Float_t factor[3]={1.,1.,1.};
 938               if(fUseStrangeness){
 939                 for(Int_t iprong=0;iprong<3;iprong++){
 940                   AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
 941                   Int_t labd= trad->GetLabel();
 942                   if(labd>=0){
 943                     AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
 944                     if(dau){
 945                       Int_t labm = dau->GetMother();
 946                       if(labm>=0){
 947                         AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
 948                         if(mot){
 949                           if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
 950                             if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
 951                             else factor[iprong]=1./.6;
 952                             //          fNentries->Fill(12);
 953                           }
 954                           if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
 955                             factor[iprong]=1./0.25;
 956                             //            fNentries->Fill(13);
 957                           }//if 3122
 958                         }//if(mot)
 959                       }//if labm>0
 960                     }//if(dau)
 961                   }//if labd>=0
 962                 }//prong loop
 963               }
 964               Float_t fact=1.;for(Int_t k=0;k<3;k++)fact=fact*factor[k];
 965               fMassHist[index]->Fill(invMass);
 966               if(fCutsDistr){
 967                 fCosPHist[index]->Fill(cosp,fact);
 968                 fDLenHist[index]->Fill(dlen,fact);
 969                 fDLxy[index]->Fill(dlxy);
 970                 fCosxy[index]->Fill(cxy);
 971
 972                 Float_t sumd02s=d->Getd0Prong(0)*d->Getd0Prong(0)*factor[0]*factor[0]+d->Getd0Prong(1)*d->Getd0Prong(1)*factor[1]*factor[1]+d->Getd0Prong(2)*d->Getd0Prong(2)*factor[2]*factor[2];
 973                 fSumd02Hist[index]->Fill(sumd02s);
 974                 fSigVertHist[index]->Fill(sigvert,fact);
 975                 fPtMaxHist[index]->Fill(ptmax,fact);
 976                 fPtKHist[index]->Fill(d->PtProng(1),fact);
 977                 fPtpi1Hist[index]->Fill(d->PtProng(0),fact);
 978                 fPtpi2Hist[index]->Fill(d->PtProng(2),fact);
 979                 fDCAHist[index]->Fill(dca,fact);
 980                 fCorreld0Kd0pi[1]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
 981                                         d->Getd0Prong(2)*d->Getd0Prong(1));
 982               }
 983               if(passTightCuts){
 984                 fMassHistTC[index]->Fill(invMass);
 985                 if(fCutsDistr){
 986                   fDLxyTC[index]->Fill(dlxy);
 987                   fCosxyTC[index]->Fill(cxy);
 988                 }
 989                 if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
 990                 else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
 991                 if(fDoImpPar){
 992                   if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse);
 993                   else{
 994                     fHistMassPtImpParTC[2]->Fill(arrayForSparse);
 995                     fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue);
 996                   }
 997                 }
 998               }
 999             }
1000             fYVsPtSig->Fill(ptCand,rapid);
1001             if(passTightCuts) fYVsPtSigTC->Fill(ptCand,rapid);
1002           }else{
1003             index=GetBackgroundHistoIndex(iPtBin);
1004             if(isFidAcc){
1005               Float_t factor[3]={1.,1.,1.};
1006               if(fUseStrangeness){
1007                 for(Int_t iprong=0;iprong<3;iprong++){
1008                   AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
1009                   Int_t labd= trad->GetLabel();
1010                   if(labd>=0){
1011                     AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
1012                     if(dau){
1013                       Int_t labm = dau->GetMother();
1014                       if(labm>=0){
1015                         AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
1016                         if(mot){
1017                           if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
1018                             if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
1019                             else factor[iprong]=1./.6;
1020                             //          fNentries->Fill(12);
1021                           }
1022                           if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
1023                             factor[iprong]=1./0.25;
1024                             //            fNentries->Fill(13);
1025                           }//if 3122
1026                         }//if(mot)
1027                       }//if labm>0
1028                     }//if(dau)
1029                   }//if labd>=0
1030                 }//prong loop
1031               }
1032
1033               Float_t fact=1.;for(Int_t k=0;k<3;k++)fact=fact*factor[k];
1034               fMassHist[index]->Fill(invMass);
1035               if(fCutsDistr){
1036                 fCosPHist[index]->Fill(cosp,fact);
1037                 fDLenHist[index]->Fill(dlen,fact);
1038                 fDLxy[index]->Fill(dlxy);
1039                 fCosxy[index]->Fill(cxy);
1040
1041                 Float_t sumd02s=d->Getd0Prong(0)*d->Getd0Prong(0)*factor[0]*factor[0]+d->Getd0Prong(1)*d->Getd0Prong(1)*factor[1]*factor[1]+d->Getd0Prong(2)*d->Getd0Prong(2)*factor[2]*factor[2];
1042                 fSumd02Hist[index]->Fill(sumd02s);
1043                 fSigVertHist[index]->Fill(sigvert,fact);
1044                 fPtMaxHist[index]->Fill(ptmax,fact);
1045                 fPtKHist[index]->Fill(d->PtProng(1),fact);
1046                 fPtpi1Hist[index]->Fill(d->PtProng(0),fact);
1047                 fPtpi2Hist[index]->Fill(d->PtProng(2),fact);
1048                 fDCAHist[index]->Fill(dca,fact);
1049                 fCorreld0Kd0pi[2]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
1050                                         d->Getd0Prong(2)*d->Getd0Prong(1));
1051               }
1052               if(passTightCuts){
1053                 fMassHistTC[index]->Fill(invMass);
1054                 if(fCutsDistr){
1055                   fDLxyTC[index]->Fill(dlxy);
1056                   fCosxyTC[index]->Fill(cxy);
1057                 }
1058                 if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
1059                 else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
1060                 if(fDoImpPar){
1061                   fHistMassPtImpParTC[4]->Fill(arrayForSparse);
1062                 }
1063               }
1064             }
1065           }
1066
1067         }
1068       }
1069
1070       if(recVtx)fRDCutsProduction->CleanOwnPrimaryVtx(d,aod,origownvtx);
1071     }
1072     if(unsetvtx) d->UnsetOwnPrimaryVtx();
1073   }
1074   fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
1075   fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
1076   }
1077   //start LS analysis
1078   if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
1079
1080   PostData(1,fOutput);
1081   PostData(2,fListCuts);
1082   PostData(3,fCounter);
1083   return;
1084 }
1085
1086 //________________________________________________________________________
1087 void AliAnalysisTaskSEDplus::CreateLikeSignHistos(){
1088   // Histos for Like Sign bckground
1089
1090   TString hisname;
1091   Int_t indexLS=0;
1092   Int_t index=0;
1093   Int_t nbins=GetNBinsHistos();
1094   for(Int_t i=0;i<fNPtBins;i++){
1095
1096     index=GetHistoIndex(i);
1097     indexLS=GetLSHistoIndex(i);
1098
1099     hisname.Form("hLSPt%dLC",i);
1100     fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1101     fMassHistLS[indexLS]->Sumw2();
1102     hisname.Form("hLSPt%dTC",i);
1103     fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1104     fMassHistLSTC[indexLS]->Sumw2();
1105
1106     hisname.Form("hCosPAllPt%dLS",i);
1107     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1108     fCosPHistLS[index]->Sumw2();
1109     hisname.Form("hDLenAllPt%dLS",i);
1110     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1111     fDLenHistLS[index]->Sumw2();
1112     hisname.Form("hSumd02AllPt%dLS",i);
1113     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1114     fSumd02HistLS[index]->Sumw2();
1115     hisname.Form("hSigVertAllPt%dLS",i);
1116     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1117     fSigVertHistLS[index]->Sumw2();
1118     hisname.Form("hPtMaxAllPt%dLS",i);
1119     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1120     fPtMaxHistLS[index]->Sumw2();
1121     hisname.Form("hDCAAllPt%dLS",i);
1122     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1123     fDCAHistLS[index]->Sumw2();
1124
1125     index=GetSignalHistoIndex(i);
1126     indexLS++;
1127
1128     hisname.Form("hLSPt%dLCnw",i);
1129     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1130     fMassHistLS[indexLS]->Sumw2();
1131     hisname.Form("hLSPt%dTCnw",i);
1132     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1133     fMassHistLSTC[indexLS]->Sumw2();
1134
1135     hisname.Form("hCosPSigPt%dLS",i);
1136     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1137     fCosPHistLS[index]->Sumw2();
1138     hisname.Form("hDLenSigPt%dLS",i);
1139     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1140     fDLenHistLS[index]->Sumw2();
1141     hisname.Form("hSumd02SigPt%dLS",i);
1142     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1143     fSumd02HistLS[index]->Sumw2();
1144     hisname.Form("hSigVertSigPt%dLS",i);
1145     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1146     fSigVertHistLS[index]->Sumw2();
1147     hisname.Form("hPtMaxSigPt%dLS",i);
1148     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1149     fPtMaxHistLS[index]->Sumw2();
1150     hisname.Form("hDCASigPt%dLS",i);
1151     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1152     fDCAHistLS[index]->Sumw2();
1153
1154     index=GetBackgroundHistoIndex(i);
1155     indexLS++;
1156
1157     hisname.Form("hLSPt%dLCntrip",i);
1158     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1159     fMassHistLS[indexLS]->Sumw2();
1160     hisname.Form("hLSPt%dTCntrip",i);
1161     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1162     fMassHistLSTC[indexLS]->Sumw2();
1163
1164     hisname.Form("hCosPBkgPt%dLS",i);
1165     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1166     fCosPHistLS[index]->Sumw2();
1167     hisname.Form("hDLenBkgPt%dLS",i);
1168     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1169     fDLenHistLS[index]->Sumw2();
1170     hisname.Form("hSumd02BkgPt%dLS",i);
1171     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1172     fSumd02HistLS[index]->Sumw2();
1173     hisname.Form("hSigVertBkgPt%dLS",i);
1174     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1175     fSigVertHistLS[index]->Sumw2();
1176     hisname.Form("hPtMaxBkgPt%dLS",i);
1177     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1178     fPtMaxHistLS[index]->Sumw2();
1179     hisname.Form("hDCABkgPt%dLS",i);
1180     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1181     fDCAHistLS[index]->Sumw2();
1182
1183     indexLS++;
1184     hisname.Form("hLSPt%dLCntripsinglecut",i);
1185     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1186     fMassHistLS[indexLS]->Sumw2();
1187     hisname.Form("hLSPt%dTCntripsinglecut",i);
1188     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1189     fMassHistLSTC[indexLS]->Sumw2();
1190
1191     indexLS++;
1192     hisname.Form("hLSPt%dLCspc",i);
1193     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1194     fMassHistLS[indexLS]->Sumw2();
1195     hisname.Form("hLSPt%dTCspc",i);
1196     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1197     fMassHistLSTC[indexLS]->Sumw2();
1198   }
1199
1200   for(Int_t i=0; i<3*fNPtBins; i++){
1201     fOutput->Add(fCosPHistLS[i]);
1202     fOutput->Add(fDLenHistLS[i]);
1203     fOutput->Add(fSumd02HistLS[i]);
1204     fOutput->Add(fSigVertHistLS[i]);
1205     fOutput->Add(fPtMaxHistLS[i]);
1206     fOutput->Add(fDCAHistLS[i]);
1207   }
1208   for(Int_t i=0; i<5*fNPtBins; i++){
1209     fOutput->Add(fMassHistLS[i]);
1210     fOutput->Add(fMassHistLSTC[i]);
1211   }
1212 }
1213
1214 //________________________________________________________________________
1215 void AliAnalysisTaskSEDplus::CreateImpactParameterHistos(){
1216   // Histos for impact paramter study
1217
1218   Int_t nmassbins=GetNBinsHistos();
1219   Int_t nbins[3]={nmassbins,200,fNImpParBins};
1220   Double_t xmin[3]={fLowmasslimit,0.,fLowerImpPar};
1221   Double_t xmax[3]={fUpmasslimit,20.,fHigherImpPar};
1222
1223   fHistMassPtImpParTC[0]=new THnSparseF("hMassPtImpParAll",
1224                                         "Mass vs. pt vs.imppar - All",
1225                                         3,nbins,xmin,xmax);
1226   fHistMassPtImpParTC[1]=new THnSparseF("hMassPtImpParPrompt",
1227                                         "Mass vs. pt vs.imppar - promptD",
1228                                         3,nbins,xmin,xmax);
1229   fHistMassPtImpParTC[2]=new THnSparseF("hMassPtImpParBfeed",
1230                                         "Mass vs. pt vs.imppar - DfromB",
1231                                         3,nbins,xmin,xmax);
1232   fHistMassPtImpParTC[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1233                                         "Mass vs. pt vs.true imppar -DfromB",
1234                                         3,nbins,xmin,xmax);
1235   fHistMassPtImpParTC[4]=new THnSparseF("hMassPtImpParBkg",
1236                                         "Mass vs. pt vs.imppar - backgr.",
1237                                         3,nbins,xmin,xmax);
1238
1239   for(Int_t i=0; i<5;i++){
1240     fOutput->Add(fHistMassPtImpParTC[i]);
1241   }
1242 }
1243
1244 //________________________________________________________________________
1245 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
1246 {
1247   // Terminate analysis
1248   //
1249   if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
1250
1251   fOutput = dynamic_cast<TList*> (GetOutputData(1));
1252   if (!fOutput) {
1253     printf("ERROR: fOutput not available\n");
1254     return;
1255   }
1256   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
1257
1258   TString hisname;
1259   Int_t index=0;
1260
1261   for(Int_t i=0;i<fNPtBins;i++){
1262     index=GetHistoIndex(i);
1263
1264     hisname.Form("hMassPt%dTC",i);
1265     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1266   }
1267
1268   TCanvas *c1=new TCanvas("c1","D+ invariant mass distribution",500,500);
1269   c1->cd();
1270   TH1F *hMassPt=(TH1F*)fOutput->FindObject("hMassPt3TC");
1271   hMassPt->SetLineColor(kBlue);
1272   hMassPt->SetXTitle("M[GeV/c^{2}]");
1273   hMassPt->Draw();
1274
1275   return;
1276 }
1277 //_________________________________________________________________________________________________
1278 Int_t AliAnalysisTaskSEDplus::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
1279   //
1280   // checking whether the mother of the particles come from a charm or a bottom quark
1281   //
1282
1283   Int_t pdgGranma = 0;
1284   Int_t mother = 0;
1285   mother = mcPartCandidate->GetMother();
1286   Int_t istep = 0;
1287   Int_t abspdgGranma =0;
1288   Bool_t isFromB=kFALSE;
1289   Bool_t isQuarkFound=kFALSE;
1290   while (mother >0 ){
1291     istep++;
1292     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1293     if (mcGranma){
1294       pdgGranma = mcGranma->GetPdgCode();
1295       abspdgGranma = TMath::Abs(pdgGranma);
1296       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1297         isFromB=kTRUE;
1298       }
1299       if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1300       mother = mcGranma->GetMother();
1301     }else{
1302       AliError("Failed casting the mother particle!");
1303       break;
1304     }
1305   }
1306
1307   if(isFromB) return 5;
1308   else return 4;
1309 }
1310 //_________________________________________________________________________________________________
1311 Float_t AliAnalysisTaskSEDplus::GetTrueImpactParameter(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partDp) const {
1312   // true impact parameter calculation
1313
1314   Double_t vtxTrue[3];
1315   mcHeader->GetVertex(vtxTrue);
1316   Double_t origD[3];
1317   partDp->XvYvZv(origD);
1318   Short_t charge=partDp->Charge();
1319   Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
1320   for(Int_t iDau=0; iDau<3; iDau++){
1321     pXdauTrue[iDau]=0.;
1322     pYdauTrue[iDau]=0.;
1323     pZdauTrue[iDau]=0.;
1324   }
1325
1326   Int_t nDau=partDp->GetNDaughters();
1327   Int_t labelFirstDau = partDp->GetDaughter(0);
1328   if(nDau==3){
1329     for(Int_t iDau=0; iDau<3; iDau++){
1330       Int_t ind = labelFirstDau+iDau;
1331       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1332       if(!part){
1333         AliError("Daughter particle not found in MC array");
1334         return 99999.;
1335       }
1336       pXdauTrue[iDau]=part->Px();
1337       pYdauTrue[iDau]=part->Py();
1338       pZdauTrue[iDau]=part->Pz();
1339     }
1340   }else if(nDau==2){
1341     Int_t theDau=0;
1342     for(Int_t iDau=0; iDau<2; iDau++){
1343       Int_t ind = labelFirstDau+iDau;
1344       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1345       if(!part){
1346         AliError("Daughter particle not found in MC array");
1347         return 99999.;
1348       }
1349       Int_t pdgCode=TMath::Abs(part->GetPdgCode());
1350       if(pdgCode==211 || pdgCode==321){
1351         pXdauTrue[theDau]=part->Px();
1352         pYdauTrue[theDau]=part->Py();
1353         pZdauTrue[theDau]=part->Pz();
1354         ++theDau;
1355       }else{
1356         Int_t nDauRes=part->GetNDaughters();
1357         if(nDauRes==2){
1358           Int_t labelFirstDauRes = part->GetDaughter(0);
1359           for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
1360             Int_t indDR = labelFirstDauRes+iDauRes;
1361             AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
1362             if(!partDR){
1363               AliError("Daughter particle not found in MC array");
1364               return 99999.;
1365             }
1366
1367             Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
1368             if(pdgCodeDR==211 || pdgCodeDR==321){
1369               pXdauTrue[theDau]=partDR->Px();
1370               pYdauTrue[theDau]=partDR->Py();
1371               pZdauTrue[theDau]=partDR->Pz();
1372               ++theDau;
1373             }
1374           }
1375         }
1376       }
1377     }
1378     if(theDau!=3){
1379       AliError("Wrong number of decay prongs");
1380       return 99999.;
1381     }
1382   }
1383
1384   Double_t d0dummy[3]={0.,0.,0.};
1385   AliAODRecoDecayHF aodDplusMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
1386   return aodDplusMC.ImpParXY();
1387
1388 }