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