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