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   fListCuts->SetOwner();
 408   AliRDHFCutsDplustoKpipi *analysis = new AliRDHFCutsDplustoKpipi(*fRDCutsAnalysis);
 409   analysis->SetName("AnalysisCuts");
 410
 411   fListCuts->Add(analysis);
 412   PostData(2,fListCuts);
 413
 414   return;
 415 }
 416
 417 //________________________________________________________________________
 418 void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
 419 {
 420   // Create the output container
 421   //
 422   if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");
 423
 424   // Several histograms are more conveniently managed in a TList
 425   fOutput = new TList();
 426   fOutput->SetOwner();
 427   fOutput->SetName("OutputHistos");
 428
 429   TString hisname;
 430   Int_t index=0;
 431   Int_t nbins=GetNBinsHistos();
 432   fHistCentrality[0]=new TH2F("hCentrMult","centrality",100,0.5,30000.5,40,0.,100.);
 433   fHistCentrality[1]=new TH2F("hCentrMult(selectedCent)","centrality(selectedCent)",100,0.5,30000.5,40,0.,100.);
 434   fHistCentrality[2]=new TH2F("hCentrMult(OutofCent)","centrality(OutofCent)",100,0.5,30000.5,40,0.,100.);
 435   for(Int_t i=0;i<3;i++){
 436     fHistCentrality[i]->Sumw2();
 437     fOutput->Add(fHistCentrality[i]);
 438   }
 439   for(Int_t i=0;i<fNPtBins;i++){
 440
 441     index=GetHistoIndex(i);
 442
 443     hisname.Form("hMassPt%d",i);
 444     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 445     fMassHist[index]->Sumw2();
 446     hisname.Form("hCosPAllPt%d",i);
 447     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
 448     fCosPHist[index]->Sumw2();
 449     hisname.Form("hDLenAllPt%d",i);
 450     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
 451     fDLenHist[index]->Sumw2();
 452     hisname.Form("hSumd02AllPt%d",i);
 453     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
 454     fSumd02Hist[index]->Sumw2();
 455     hisname.Form("hSigVertAllPt%d",i);
 456     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 457     fSigVertHist[index]->Sumw2();
 458     hisname.Form("hPtMaxAllPt%d",i);
 459     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
 460     fPtMaxHist[index]->Sumw2();
 461     hisname.Form("hPtKPt%d",i);
 462     fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 463     fPtKHist[index]->Sumw2();
 464     hisname.Form("hPtpi1Pt%d",i);
 465     fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 466     fPtpi1Hist[index]->Sumw2();
 467     hisname.Form("hPtpi2Pt%d",i);
 468     fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 469     fPtpi2Hist[index]->Sumw2();
 470     hisname.Form("hDCAAllPt%d",i);
 471     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 472     fDCAHist[index]->Sumw2();
 473
 474     hisname.Form("hDLxyPt%d",i);
 475     fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
 476     fDLxy[index]->Sumw2();
 477     hisname.Form("hCosxyPt%d",i);
 478     fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
 479     fCosxy[index]->Sumw2();
 480     hisname.Form("hDLxyPt%dTC",i);
 481     fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
 482     fDLxyTC[index]->Sumw2();
 483     hisname.Form("hCosxyPt%dTC",i);
 484     fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
 485     fCosxyTC[index]->Sumw2();
 486
 487     hisname.Form("hMassPt%dTC",i);
 488     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 489     fMassHistTC[index]->Sumw2();
 490     hisname.Form("hMassPt%dTCPlus",i);
 491     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 492     fMassHistTCPlus[index]->Sumw2();
 493     hisname.Form("hMassPt%dTCMinus",i);
 494     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 495     fMassHistTCMinus[index]->Sumw2();
 496
 497
 498
 499     index=GetSignalHistoIndex(i);
 500     hisname.Form("hSigPt%d",i);
 501     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 502     fMassHist[index]->Sumw2();
 503     hisname.Form("hCosPSigPt%d",i);
 504     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
 505     fCosPHist[index]->Sumw2();
 506     hisname.Form("hDLenSigPt%d",i);
 507     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
 508     fDLenHist[index]->Sumw2();
 509     hisname.Form("hSumd02SigPt%d",i);
 510     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
 511     fSumd02Hist[index]->Sumw2();
 512     hisname.Form("hSigVertSigPt%d",i);
 513     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 514     fSigVertHist[index]->Sumw2();
 515     hisname.Form("hPtMaxSigPt%d",i);
 516     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
 517     fPtMaxHist[index]->Sumw2();
 518     hisname.Form("hPtKSigPt%d",i);
 519     fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 520     fPtKHist[index]->Sumw2();
 521     hisname.Form("hPtpi1SigPt%d",i);
 522     fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 523     fPtpi1Hist[index]->Sumw2();
 524     hisname.Form("hPtpi2SigPt%d",i);
 525     fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 526     fPtpi2Hist[index]->Sumw2();
 527
 528     hisname.Form("hDCASigPt%d",i);
 529     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 530     fDCAHist[index]->Sumw2();
 531
 532     hisname.Form("hDLxySigPt%d",i);
 533     fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
 534     fDLxy[index]->Sumw2();
 535     hisname.Form("hCosxySigPt%d",i);
 536     fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
 537     fCosxy[index]->Sumw2();
 538     hisname.Form("hDLxySigPt%dTC",i);
 539     fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
 540     fDLxyTC[index]->Sumw2();
 541     hisname.Form("hCosxySigPt%dTC",i);
 542     fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
 543     fCosxyTC[index]->Sumw2();
 544     hisname.Form("hSigPt%dTC",i);
 545     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 546     fMassHistTC[index]->Sumw2();
 547     hisname.Form("hSigPt%dTCPlus",i);
 548     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 549     fMassHistTCPlus[index]->Sumw2();
 550     hisname.Form("hSigPt%dTCMinus",i);
 551     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 552     fMassHistTCMinus[index]->Sumw2();
 553
 554
 555     index=GetBackgroundHistoIndex(i);
 556     hisname.Form("hBkgPt%d",i);
 557     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 558     fMassHist[index]->Sumw2();
 559     hisname.Form("hCosPBkgPt%d",i);
 560     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,1.);
 561     fCosPHist[index]->Sumw2();
 562     hisname.Form("hDLenBkgPt%d",i);
 563     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.5);
 564     fDLenHist[index]->Sumw2();
 565     hisname.Form("hSumd02BkgPt%d",i);
 566     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,1.);
 567     fSumd02Hist[index]->Sumw2();
 568     hisname.Form("hSigVertBkgPt%d",i);
 569     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 570     fSigVertHist[index]->Sumw2();
 571     hisname.Form("hPtMaxBkgPt%d",i);
 572     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.5,5.);
 573     fPtMaxHist[index]->Sumw2();
 574     hisname.Form("hPtKBkgPt%d",i);
 575     fPtKHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 576     fPtKHist[index]->Sumw2();
 577     hisname.Form("hPtpi1BkgPt%d",i);
 578     fPtpi1Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 579     fPtpi1Hist[index]->Sumw2();
 580     hisname.Form("hPtpi2BkgPt%d",i);
 581     fPtpi2Hist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,5.);
 582     fPtpi2Hist[index]->Sumw2();
 583     hisname.Form("hDCABkgPt%d",i);
 584     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,0.1);
 585     fDCAHist[index]->Sumw2();
 586
 587     hisname.Form("hDLxyBkgPt%d",i);
 588     fDLxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
 589     fDLxy[index]->Sumw2();
 590     hisname.Form("hCosxyBkgPt%d",i);
 591     fCosxy[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
 592     fCosxy[index]->Sumw2();
 593     hisname.Form("hDLxyBkgPt%dTC",i);
 594     fDLxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,0.,10.);
 595     fDLxyTC[index]->Sumw2();
 596     hisname.Form("hCosxyBkgPt%dTC",i);
 597     fCosxyTC[index]=new TH1F(hisname.Data(),hisname.Data(),100,-1,1.);
 598     fCosxyTC[index]->Sumw2();
 599
 600
 601     hisname.Form("hBkgPt%dTC",i);
 602     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 603     fMassHistTC[index]->Sumw2();
 604     hisname.Form("hBkgPt%dTCPlus",i);
 605     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 606     fMassHistTCPlus[index]->Sumw2();
 607     hisname.Form("hBkgPt%dTCMinus",i);
 608     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 609     fMassHistTCMinus[index]->Sumw2();
 610   }
 611
 612
 613   for(Int_t i=0; i<3*fNPtBins; i++){
 614     fOutput->Add(fMassHist[i]);
 615     if(fCutsDistr){
 616       fOutput->Add(fCosPHist[i]);
 617       fOutput->Add(fDLenHist[i]);
 618       fOutput->Add(fSumd02Hist[i]);
 619       fOutput->Add(fSigVertHist[i]);
 620       fOutput->Add(fPtMaxHist[i]);
 621       fOutput->Add(fPtKHist[i]);
 622       fOutput->Add(fPtpi1Hist[i]);
 623       fOutput->Add(fPtpi2Hist[i]);
 624       fOutput->Add(fDCAHist[i]);
 625       fOutput->Add(fDLxy[i]);
 626       fOutput->Add(fDLxyTC[i]);
 627       fOutput->Add(fCosxy[i]);
 628       fOutput->Add(fCosxyTC[i]);
 629     }
 630     fOutput->Add(fMassHistTC[i]);
 631     fOutput->Add(fMassHistTCPlus[i]);
 632     fOutput->Add(fMassHistTCMinus[i]);
 633
 634   }
 635
 636   if(fCutsDistr){
 637     fCorreld0Kd0pi[0]=new TH2F("hCorreld0Kd0piAll","",100,-0.02,0.02,100,-0.02,0.02);
 638     fCorreld0Kd0pi[1]=new TH2F("hCorreld0Kd0piSig","",100,-0.02,0.02,100,-0.02,0.02);
 639     fCorreld0Kd0pi[2]=new TH2F("hCorreld0Kd0piBkg","",100,-0.02,0.02,100,-0.02,0.02);
 640     for(Int_t i=0; i<3; i++){
 641       fCorreld0Kd0pi[i]->Sumw2();
 642       fOutput->Add(fCorreld0Kd0pi[i]);
 643     }
 644   }
 645
 646   fHistNEvents = new TH1F("fHistNEvents", "number of events ",11,-0.5,10.5);
 647   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 648   fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents accepted");
 649   fHistNEvents->GetXaxis()->SetBinLabel(3,"Rejected due to trigger");
 650   fHistNEvents->GetXaxis()->SetBinLabel(4,"Rejected pileup events");
 651   fHistNEvents->GetXaxis()->SetBinLabel(5,"Rejected due to centrality");
 652   fHistNEvents->GetXaxis()->SetBinLabel(6,"Rejected due to vtxz");
 653   fHistNEvents->GetXaxis()->SetBinLabel(7,"Rejected due to Physics Sel");
 654   fHistNEvents->GetXaxis()->SetBinLabel(8,"Total no. of candidate");
 655   fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
 656   fHistNEvents->GetXaxis()->SetBinLabel(10,"D+ after loose cuts");
 657   fHistNEvents->GetXaxis()->SetBinLabel(11,"D+ after tight cuts");
 658
 659   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
 660   fHistNEvents->Sumw2();
 661   fHistNEvents->SetMinimum(0);
 662   fOutput->Add(fHistNEvents);
 663
 664   fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
 665   fPtVsMassTC=new TH2F("hPtVsMassTC","PtVsMass (analysis cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
 666   fYVsPt=new TH3F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.,nbins,fLowmasslimit,fUpmasslimit);
 667   fYVsPtTC=new TH3F("hYVsPtTC","YvsPt (analysis cuts)",40,0.,20.,80,-2.,2.,nbins,fLowmasslimit,fUpmasslimit);
 668   fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
 669   fYVsPtSigTC=new TH2F("hYVsPtSigTC","YvsPt (MC, only Sig, analysis cuts)",40,0.,20.,80,-2.,2.);
 670   fPhiEtaCand=new TH2F("hPhiEtaCand","phi vs. eta candidates",20,-1.,1.,50,0.,2*TMath::Pi());
 671   fPhiEtaCandSigReg=new TH2F("hPhiEtaCandSigReg","phi vs. eta candidates",20,-1.,1.,50,0.,2*TMath::Pi());
 672   fSPDMult = new TH1F("hSPDMult", "Tracklets multiplicity; Tracklets ; Entries",200,0.,200.);
 673   fOutput->Add(fPtVsMass);
 674   fOutput->Add(fPtVsMassTC);
 675   fOutput->Add(fYVsPt);
 676   fOutput->Add(fYVsPtTC);
 677   fOutput->Add(fYVsPtSig);
 678   fOutput->Add(fYVsPtSigTC);
 679   fOutput->Add(fPhiEtaCand);
 680   fOutput->Add(fPhiEtaCandSigReg);
 681   fOutput->Add(fSPDMult);
 682
 683
 684   //Counter for Normalization
 685   TString normName="NormalizationCounter";
 686   AliAnalysisDataContainer *cont = GetOutputSlot(3)->GetContainer();
 687   if(cont)normName=(TString)cont->GetName();
 688   fCounter = new AliNormalizationCounter(normName.Data());
 689   fCounter->Init();
 690
 691   if(fDoLS) CreateLikeSignHistos();
 692   if(fDoImpPar) CreateImpactParameterHistos();
 693
 694   if(fFillNtuple){
 695     OpenFile(4); // 4 is the slot number of the ntuple
 696
 697
 698     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");
 699
 700   }
 701
 702   return;
 703 }
 704
 705 //________________________________________________________________________
 706 void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
 707 {
 708   // Execute analysis for current event:
 709   // heavy flavor candidates association to MC truth
 710
 711   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 712
 713   TClonesArray *array3Prong = 0;
 714   TClonesArray *arrayLikeSign =0;
 715   if(!aod && AODEvent() && IsStandardAOD()) {
 716     // In case there is an AOD handler writing a standard AOD, use the AOD
 717     // event in memory rather than the input (ESD) event.
 718     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 719     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 720     // have to taken from the AOD event hold by the AliAODExtension
 721     AliAODHandler* aodHandler = (AliAODHandler*)
 722       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 723     if(aodHandler->GetExtensions()) {
 724       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 725       AliAODEvent *aodFromExt = ext->GetAOD();
 726       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 727       arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
 728     }
 729   } else if(aod) {
 730     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 731     arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
 732   }
 733
 734   if(!aod || !array3Prong) {
 735     printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
 736     return;
 737   }
 738   if(!arrayLikeSign && fDoLS) {
 739     printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
 740     return;
 741   }
 742
 743   // fix for temporary bug in ESDfilter
 744   // the AODs with null vertex pointer didn't pass the PhysSel
 745   if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
 746
 747   //Store the event in AliNormalizationCounter->To disable for Pb-Pb? Add a flag to disable it?
 748   fCounter->StoreEvent(aod,fRDCutsAnalysis,fReadMC);
 749
 750   fHistNEvents->Fill(0); // count event
 751   Int_t runNumber=aod->GetRunNumber();
 752
 753   //Event selection
 754   Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
 755   Float_t ntracks=aod->GetNTracks();
 756   Float_t evCentr=0;
 757   if(fRDCutsAnalysis->GetUseCentrality()>0) evCentr=fRDCutsAnalysis->GetCentrality(aod);
 758   fHistCentrality[0]->Fill(ntracks,evCentr);
 759   if(fRDCutsAnalysis->GetWhyRejection()==5) fHistNEvents->Fill(2);
 760   if(fRDCutsAnalysis->GetWhyRejection()==1) fHistNEvents->Fill(3);
 761   if(fRDCutsAnalysis->GetWhyRejection()==2){fHistNEvents->Fill(4);fHistCentrality[2]->Fill(ntracks,evCentr);}
 762   if(fRDCutsAnalysis->GetWhyRejection()==6)fHistNEvents->Fill(5);
 763   if(fRDCutsAnalysis->GetWhyRejection()==7)fHistNEvents->Fill(6);
 764
 765   // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
 766   //TString trigclass=aod->GetFiredTriggerClasses();
 767   // Post the data already here
 768   PostData(1,fOutput);
 769   if(!isEvSel)return;
 770   Int_t tracklets=AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aod,-1.,1.);
 771   // printf("ntracklet===%d\n",tracklets);
 772   fSPDMult->Fill(tracklets);
 773
 774   fHistCentrality[1]->Fill(ntracks,evCentr);
 775   fHistNEvents->Fill(1);
 776
 777   TClonesArray *arrayMC=0;
 778   AliAODMCHeader *mcHeader=0;
 779
 780   // AOD primary vertex
 781   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 782   //    vtx1->Print();
 783   //   TString primTitle = vtx1->GetTitle();
 784   //if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0)fHistNEvents->Fill(2);
 785
 786   // load MC particles
 787   if(fReadMC){
 788     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 789     if(!arrayMC) {
 790       printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
 791       return;
 792     }
 793
 794     // load MC header
 795     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 796     if(!mcHeader) {
 797       printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
 798       return;
 799     }
 800   }
 801
 802   Int_t n3Prong = array3Prong->GetEntriesFast();
 803   //  printf("Number of D+->Kpipi: %d and of tracks: %d\n",n3Prong,aod->GetNumberOfTracks());
 804
 805   Int_t nOS=0;
 806   Int_t index;
 807   Int_t pdgDgDplustoKpipi[3]={321,211,211};
 808
 809   if(fDoLS>1){//Normalizations for LS
 810     for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 811       AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 812       if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
 813         if(fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod))nOS++;
 814       }
 815     }
 816   }else{//Standard analysis
 817     // 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
 818     //Double_t *cutsDplus = new (Double_t*)fRDCuts->GetCuts();
 819     Int_t nSelectedloose=0,nSelectedtight=0;
 820     for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 821       AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 822       fHistNEvents->Fill(7);
 823       if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
 824         fHistNEvents->Fill(8);
 825         continue;
 826       }
 827
 828       Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kAll,aod);
 829
 830       if(!fRDCutsAnalysis->GetIsSelectedCuts()) continue;
 831
 832       Double_t etaD=d->Eta();
 833       Double_t phiD=d->Phi();
 834       if(fEtaSelection!=0){
 835         if(fEtaSelection==1 && etaD<0) continue;
 836         if(fEtaSelection==-1 && etaD>0) continue;
 837       }
 838
 839       Bool_t unsetvtx=kFALSE;
 840       if(!d->GetOwnPrimaryVtx()){
 841         d->SetOwnPrimaryVtx(vtx1);
 842         unsetvtx=kTRUE;
 843       }
 844
 845       Double_t ptCand = d->Pt();
 846       Int_t iPtBin = fRDCutsAnalysis->PtBin(ptCand);
 847
 848       Bool_t recVtx=kFALSE;
 849       AliAODVertex *origownvtx=0x0;
 850       if(fRDCutsAnalysis->GetIsPrimaryWithoutDaughters()){
 851         if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
 852         if(fRDCutsAnalysis->RecalcOwnPrimaryVtx(d,aod))recVtx=kTRUE;
 853         else fRDCutsAnalysis->CleanOwnPrimaryVtx(d,aod,origownvtx);
 854       }
 855
 856       Int_t labDp=-1;
 857       Bool_t isPrimary=kTRUE;
 858       Float_t pdgCode=-2;
 859       Float_t trueImpParXY=0.;
 860       if(fReadMC){
 861         labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
 862         if(labDp>=0){
 863           AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
 864           if(CheckOrigin(arrayMC,partDp)==5) isPrimary=kFALSE;
 865           pdgCode=TMath::Abs(partDp->GetPdgCode());
 866           if(!isPrimary){
 867             trueImpParXY=GetTrueImpactParameter(mcHeader,arrayMC,partDp)*10000.;
 868           }
 869         }else{
 870           pdgCode=-1;
 871         }
 872       }
 873
 874       Double_t invMass=d->InvMassDplus();
 875       Double_t rapid=d->YDplus();
 876       fYVsPt->Fill(ptCand,rapid,invMass);
 877       if(passTightCuts) {fYVsPtTC->Fill(ptCand,rapid,invMass);nOS++;}
 878       Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
 879       if(isFidAcc){
 880         fPtVsMass->Fill(invMass,ptCand);
 881         if(passTightCuts){
 882           fPtVsMassTC->Fill(invMass,ptCand);
 883           fPhiEtaCand->Fill(etaD,phiD);
 884           if(TMath::Abs(invMass-1.8696)<0.05) fPhiEtaCandSigReg->Fill(etaD,phiD);
 885         }
 886       }
 887
 888
 889       Double_t  dlen=0,cosp=0,maxdca=0,sigvert=0,sumD02=0,ptmax=0,dlxy=0,cxy=0;
 890       if(fCutsDistr||fFillNtuple||fDoImpPar){
 891         dlen=d->DecayLength();
 892         cosp=d->CosPointingAngle();
 893         sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 894         maxdca=-9999.;
 895         for(Int_t idau=0;idau<3;idau++) if(d->GetDCA(idau)>maxdca) maxdca=d->GetDCA(idau);
 896         sigvert=d->GetSigmaVert();
 897         ptmax=0;
 898         for(Int_t i=0;i<3;i++){
 899           if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
 900         }
 901         dlxy=d->NormalizedDecayLengthXY();
 902         cxy=d->CosPointingAngleXY();
 903       }
 904       Double_t impparXY=d->ImpParXY()*10000.;
 905       Double_t arrayForSparse[6]={invMass,ptCand,impparXY,cosp,dlen,static_cast<Double_t>(tracklets)};
 906       Double_t arrayForSparseTrue[6]={invMass,ptCand,trueImpParXY,cosp,dlen,static_cast<Double_t>(tracklets)};
 907
 908       //Ntuple
 909       Float_t tmp[31];
 910       if(fFillNtuple){
 911         tmp[0]=pdgCode;
 912         if(!isPrimary) tmp[0]+=5000.;
 913         tmp[1]=d->Px();
 914         tmp[2]=d->Py();
 915         tmp[3]=d->Pz();
 916         tmp[4]=d->Pt();
 917         tmp[5]=d->GetCharge();
 918         //      tmp[5]=fRDCutsAnalysis->GetPIDBitMask(d);
 919         tmp[6]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(0));
 920         tmp[7]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(1));
 921         tmp[8]=fRDCutsAnalysis->GetPIDTrackTPCTOFBitMap((AliAODTrack*)d->GetDaughter(2));
 922         tmp[9]=d->PtProng(0);
 923         tmp[10]=d->PtProng(1);
 924         tmp[11]=d->PtProng(2);
 925         tmp[12]=d->PProng(0);
 926         tmp[13]=d->PProng(1);
 927         tmp[14]=d->PProng(2);
 928         tmp[15]=cosp;
 929         tmp[16]=cxy;
 930         tmp[17]=dlen;
 931         tmp[18]=d->NormalizedDecayLength();
 932         tmp[19]=d->DecayLengthXY();
 933         tmp[20]=dlxy;
 934         tmp[21]=d->InvMassDplus();
 935         tmp[22]=sigvert;
 936         tmp[23]=d->Getd0Prong(0);
 937         tmp[24]=d->Getd0Prong(1);
 938         tmp[25]=d->Getd0Prong(2);
 939         tmp[26]=maxdca;
 940         tmp[27]=ntracks;
 941         tmp[28]=fRDCutsAnalysis->GetCentrality(aod);
 942         tmp[29]=runNumber;
 943         tmp[30]=d->HasBadDaughters();
 944         fNtupleDplus->Fill(tmp);
 945         PostData(4,fNtupleDplus);
 946       }
 947
 948       //Fill histos
 949       index=GetHistoIndex(iPtBin);
 950       if(isFidAcc){
 951         fHistNEvents->Fill(9);
 952         nSelectedloose++;
 953         fMassHist[index]->Fill(invMass);
 954         if(fCutsDistr){
 955           fCosPHist[index]->Fill(cosp);
 956           fDLenHist[index]->Fill(dlen);
 957           fSumd02Hist[index]->Fill(sumD02);
 958           fSigVertHist[index]->Fill(sigvert);
 959           fPtMaxHist[index]->Fill(ptmax);
 960           fPtKHist[index]->Fill(d->PtProng(1));
 961           fPtpi1Hist[index]->Fill(d->PtProng(0));
 962           fPtpi2Hist[index]->Fill(d->PtProng(2));
 963           fDCAHist[index]->Fill(maxdca);
 964           fDLxy[index]->Fill(dlxy);
 965           fCosxy[index]->Fill(cxy);
 966           fCorreld0Kd0pi[0]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
 967                                   d->Getd0Prong(2)*d->Getd0Prong(1));
 968         }
 969         if(passTightCuts){
 970           fHistNEvents->Fill(10);
 971           nSelectedtight++;
 972           fMassHistTC[index]->Fill(invMass);
 973           if(fCutsDistr){
 974             fDLxyTC[index]->Fill(dlxy);
 975             fCosxyTC[index]->Fill(cxy);
 976           }
 977           if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
 978           else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
 979           if(fDoImpPar){
 980             fHistMassPtImpParTC[0]->Fill(arrayForSparse);
 981           }
 982         }
 983       }
 984
 985       if(fReadMC){
 986         if(isFidAcc){
 987           if(labDp>=0) {
 988             index=GetSignalHistoIndex(iPtBin);
 989             if(passTightCuts&&fDoImpPar){
 990               if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse);
 991               else{
 992                 fHistMassPtImpParTC[2]->Fill(arrayForSparse);
 993                 fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue);
 994               }
 995             }
 996           }else{
 997             index=GetBackgroundHistoIndex(iPtBin);
 998             if(passTightCuts&&fDoImpPar)fHistMassPtImpParTC[4]->Fill(arrayForSparse);
 999           }
1000
1001           fMassHist[index]->Fill(invMass);
1002           if(fCutsDistr){
1003             Float_t fact=1.;
1004             Float_t factor[3]={1.,1.,1.};
1005             if(fUseStrangeness) fact=GetStrangenessWeights(d,arrayMC,factor);
1006             fCosPHist[index]->Fill(cosp,fact);
1007             fDLenHist[index]->Fill(dlen,fact);
1008             fDLxy[index]->Fill(dlxy);
1009             fCosxy[index]->Fill(cxy);
1010
1011             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];
1012             fSumd02Hist[index]->Fill(sumd02s);
1013             fSigVertHist[index]->Fill(sigvert,fact);
1014             fPtMaxHist[index]->Fill(ptmax,fact);
1015             fPtKHist[index]->Fill(d->PtProng(1),fact);
1016             fPtpi1Hist[index]->Fill(d->PtProng(0),fact);
1017             fPtpi2Hist[index]->Fill(d->PtProng(2),fact);
1018             fDCAHist[index]->Fill(maxdca,fact);
1019             fCorreld0Kd0pi[1]->Fill(d->Getd0Prong(0)*d->Getd0Prong(1),
1020                                     d->Getd0Prong(2)*d->Getd0Prong(1));
1021           }
1022           if(passTightCuts){
1023             fMassHistTC[index]->Fill(invMass);
1024             if(fCutsDistr){
1025               fDLxyTC[index]->Fill(dlxy);
1026               fCosxyTC[index]->Fill(cxy);
1027             }
1028             if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
1029             else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
1030           }
1031         }else{//outside fidAcc
1032           if(labDp>=0){
1033             fYVsPtSig->Fill(ptCand,rapid, invMass);
1034             if(passTightCuts)fYVsPtSigTC->Fill(ptCand,rapid, invMass);
1035           }
1036         }
1037       }//readmc
1038
1039       if(recVtx)fRDCutsAnalysis->CleanOwnPrimaryVtx(d,aod,origownvtx);
1040
1041       if(unsetvtx) d->UnsetOwnPrimaryVtx();
1042     }
1043     fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
1044     fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
1045   }
1046   //start LS analysis
1047   if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
1048
1049   PostData(1,fOutput);
1050   PostData(2,fListCuts);
1051   PostData(3,fCounter);
1052   return;
1053 }
1054
1055 //________________________________________________________________________
1056 void AliAnalysisTaskSEDplus::CreateLikeSignHistos(){
1057   // Histos for Like Sign bckground
1058
1059   TString hisname;
1060   Int_t indexLS=0;
1061   Int_t index=0;
1062   Int_t nbins=GetNBinsHistos();
1063   for(Int_t i=0;i<fNPtBins;i++){
1064
1065     index=GetHistoIndex(i);
1066     indexLS=GetLSHistoIndex(i);
1067
1068     hisname.Form("hLSPt%dLC",i);
1069     fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1070     fMassHistLS[indexLS]->Sumw2();
1071     hisname.Form("hLSPt%dTC",i);
1072     fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1073     fMassHistLSTC[indexLS]->Sumw2();
1074
1075     hisname.Form("hCosPAllPt%dLS",i);
1076     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1077     fCosPHistLS[index]->Sumw2();
1078     hisname.Form("hDLenAllPt%dLS",i);
1079     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1080     fDLenHistLS[index]->Sumw2();
1081     hisname.Form("hSumd02AllPt%dLS",i);
1082     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1083     fSumd02HistLS[index]->Sumw2();
1084     hisname.Form("hSigVertAllPt%dLS",i);
1085     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1086     fSigVertHistLS[index]->Sumw2();
1087     hisname.Form("hPtMaxAllPt%dLS",i);
1088     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1089     fPtMaxHistLS[index]->Sumw2();
1090     hisname.Form("hDCAAllPt%dLS",i);
1091     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1092     fDCAHistLS[index]->Sumw2();
1093
1094     index=GetSignalHistoIndex(i);
1095     indexLS++;
1096
1097     hisname.Form("hLSPt%dLCnw",i);
1098     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1099     fMassHistLS[indexLS]->Sumw2();
1100     hisname.Form("hLSPt%dTCnw",i);
1101     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1102     fMassHistLSTC[indexLS]->Sumw2();
1103
1104     hisname.Form("hCosPSigPt%dLS",i);
1105     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1106     fCosPHistLS[index]->Sumw2();
1107     hisname.Form("hDLenSigPt%dLS",i);
1108     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1109     fDLenHistLS[index]->Sumw2();
1110     hisname.Form("hSumd02SigPt%dLS",i);
1111     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1112     fSumd02HistLS[index]->Sumw2();
1113     hisname.Form("hSigVertSigPt%dLS",i);
1114     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1115     fSigVertHistLS[index]->Sumw2();
1116     hisname.Form("hPtMaxSigPt%dLS",i);
1117     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1118     fPtMaxHistLS[index]->Sumw2();
1119     hisname.Form("hDCASigPt%dLS",i);
1120     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1121     fDCAHistLS[index]->Sumw2();
1122
1123     index=GetBackgroundHistoIndex(i);
1124     indexLS++;
1125
1126     hisname.Form("hLSPt%dLCntrip",i);
1127     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1128     fMassHistLS[indexLS]->Sumw2();
1129     hisname.Form("hLSPt%dTCntrip",i);
1130     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1131     fMassHistLSTC[indexLS]->Sumw2();
1132
1133     hisname.Form("hCosPBkgPt%dLS",i);
1134     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
1135     fCosPHistLS[index]->Sumw2();
1136     hisname.Form("hDLenBkgPt%dLS",i);
1137     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
1138     fDLenHistLS[index]->Sumw2();
1139     hisname.Form("hSumd02BkgPt%dLS",i);
1140     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
1141     fSumd02HistLS[index]->Sumw2();
1142     hisname.Form("hSigVertBkgPt%dLS",i);
1143     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1144     fSigVertHistLS[index]->Sumw2();
1145     hisname.Form("hPtMaxBkgPt%dLS",i);
1146     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
1147     fPtMaxHistLS[index]->Sumw2();
1148     hisname.Form("hDCABkgPt%dLS",i);
1149     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
1150     fDCAHistLS[index]->Sumw2();
1151
1152     indexLS++;
1153     hisname.Form("hLSPt%dLCntripsinglecut",i);
1154     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1155     fMassHistLS[indexLS]->Sumw2();
1156     hisname.Form("hLSPt%dTCntripsinglecut",i);
1157     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1158     fMassHistLSTC[indexLS]->Sumw2();
1159
1160     indexLS++;
1161     hisname.Form("hLSPt%dLCspc",i);
1162     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1163     fMassHistLS[indexLS]->Sumw2();
1164     hisname.Form("hLSPt%dTCspc",i);
1165     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
1166     fMassHistLSTC[indexLS]->Sumw2();
1167   }
1168
1169   for(Int_t i=0; i<3*fNPtBins; i++){
1170     fOutput->Add(fCosPHistLS[i]);
1171     fOutput->Add(fDLenHistLS[i]);
1172     fOutput->Add(fSumd02HistLS[i]);
1173     fOutput->Add(fSigVertHistLS[i]);
1174     fOutput->Add(fPtMaxHistLS[i]);
1175     fOutput->Add(fDCAHistLS[i]);
1176   }
1177   for(Int_t i=0; i<5*fNPtBins; i++){
1178     fOutput->Add(fMassHistLS[i]);
1179     fOutput->Add(fMassHistLSTC[i]);
1180   }
1181 }
1182
1183 //________________________________________________________________________
1184 void AliAnalysisTaskSEDplus::CreateImpactParameterHistos(){
1185   // Histos for impact paramter study
1186
1187   Int_t nmassbins=GetNBinsHistos();
1188   Double_t maxmult;
1189   if(fSystem==1) maxmult=5000;
1190   else maxmult=200;
1191   Int_t nbins[6]={nmassbins,200,fNImpParBins,5,50,100};
1192   Double_t xmin[6]={fLowmasslimit,0.,fLowerImpPar,0.95,0.,-0.5};
1193   Double_t xmax[6]={fUpmasslimit,40.,fHigherImpPar,1.,1.,maxmult};
1194
1195   fHistMassPtImpParTC[0]=new THnSparseF("hMassPtImpParAll",
1196                                         "Mass vs. pt vs.imppar - All",
1197                                         6,nbins,xmin,xmax);
1198   fHistMassPtImpParTC[1]=new THnSparseF("hMassPtImpParPrompt",
1199                                         "Mass vs. pt vs.imppar - promptD",
1200                                         6,nbins,xmin,xmax);
1201   fHistMassPtImpParTC[2]=new THnSparseF("hMassPtImpParBfeed",
1202                                         "Mass vs. pt vs.imppar - DfromB",
1203                                         6,nbins,xmin,xmax);
1204   fHistMassPtImpParTC[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1205                                         "Mass vs. pt vs.true imppar -DfromB",
1206                                         6,nbins,xmin,xmax);
1207   fHistMassPtImpParTC[4]=new THnSparseF("hMassPtImpParBkg",
1208                                         "Mass vs. pt vs.imppar - backgr.",
1209                                         6,nbins,xmin,xmax);
1210   for(Int_t i=0; i<5;i++){
1211     fOutput->Add(fHistMassPtImpParTC[i]);
1212   }
1213 }
1214
1215 //________________________________________________________________________
1216 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
1217 {
1218   // Terminate analysis
1219   //
1220   if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
1221
1222   fOutput = dynamic_cast<TList*> (GetOutputData(1));
1223   if (!fOutput) {
1224     printf("ERROR: fOutput not available\n");
1225     return;
1226   }
1227
1228   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
1229   if(fHistNEvents){
1230     printf("Number of analyzed events = %d\n",(Int_t)fHistNEvents->GetBinContent(2));
1231   }else{
1232     printf("ERROR: fHistNEvents not available\n");
1233     return;
1234   }
1235
1236   return;
1237 }
1238 //_________________________________________________________________________________________________
1239 Int_t AliAnalysisTaskSEDplus::CheckOrigin(TClonesArray* arrayMC, const AliAODMCParticle *mcPartCandidate) const {
1240   //
1241   // checking whether the mother of the particles come from a charm or a bottom quark
1242   //
1243
1244   Int_t pdgGranma = 0;
1245   Int_t mother = 0;
1246   mother = mcPartCandidate->GetMother();
1247   Int_t istep = 0;
1248   Int_t abspdgGranma =0;
1249   Bool_t isFromB=kFALSE;
1250   //  Bool_t isQuarkFound=kFALSE;
1251   while (mother >0 ){
1252     istep++;
1253     AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1254     if (mcGranma){
1255       pdgGranma = mcGranma->GetPdgCode();
1256       abspdgGranma = TMath::Abs(pdgGranma);
1257       if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1258         isFromB=kTRUE;
1259       }
1260       //     if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
1261       mother = mcGranma->GetMother();
1262     }else{
1263       AliError("Failed casting the mother particle!");
1264       break;
1265     }
1266   }
1267
1268   if(isFromB) return 5;
1269   else return 4;
1270 }
1271 //_________________________________________________________________________________________________
1272 Float_t AliAnalysisTaskSEDplus::GetTrueImpactParameter(const AliAODMCHeader *mcHeader, TClonesArray* arrayMC, const AliAODMCParticle *partDp) const {
1273   // true impact parameter calculation
1274
1275   Double_t vtxTrue[3];
1276   mcHeader->GetVertex(vtxTrue);
1277   Double_t origD[3];
1278   partDp->XvYvZv(origD);
1279   Short_t charge=partDp->Charge();
1280   Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
1281   for(Int_t iDau=0; iDau<3; iDau++){
1282     pXdauTrue[iDau]=0.;
1283     pYdauTrue[iDau]=0.;
1284     pZdauTrue[iDau]=0.;
1285   }
1286
1287   Int_t nDau=partDp->GetNDaughters();
1288   Int_t labelFirstDau = partDp->GetDaughter(0);
1289   if(nDau==3){
1290     for(Int_t iDau=0; iDau<3; iDau++){
1291       Int_t ind = labelFirstDau+iDau;
1292       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1293       if(!part){
1294         AliError("Daughter particle not found in MC array");
1295         return 99999.;
1296       }
1297       pXdauTrue[iDau]=part->Px();
1298       pYdauTrue[iDau]=part->Py();
1299       pZdauTrue[iDau]=part->Pz();
1300     }
1301   }else if(nDau==2){
1302     Int_t theDau=0;
1303     for(Int_t iDau=0; iDau<2; iDau++){
1304       Int_t ind = labelFirstDau+iDau;
1305       AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1306       if(!part){
1307         AliError("Daughter particle not found in MC array");
1308         return 99999.;
1309       }
1310       Int_t pdgCode=TMath::Abs(part->GetPdgCode());
1311       if(pdgCode==211 || pdgCode==321){
1312         pXdauTrue[theDau]=part->Px();
1313         pYdauTrue[theDau]=part->Py();
1314         pZdauTrue[theDau]=part->Pz();
1315         ++theDau;
1316       }else{
1317         Int_t nDauRes=part->GetNDaughters();
1318         if(nDauRes==2){
1319           Int_t labelFirstDauRes = part->GetDaughter(0);
1320           for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
1321             Int_t indDR = labelFirstDauRes+iDauRes;
1322             AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
1323             if(!partDR){
1324               AliError("Daughter particle not found in MC array");
1325               return 99999.;
1326             }
1327
1328             Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
1329             if(pdgCodeDR==211 || pdgCodeDR==321){
1330               pXdauTrue[theDau]=partDR->Px();
1331               pYdauTrue[theDau]=partDR->Py();
1332               pZdauTrue[theDau]=partDR->Pz();
1333               ++theDau;
1334             }
1335           }
1336         }
1337       }
1338     }
1339     if(theDau!=3){
1340       AliError("Wrong number of decay prongs");
1341       return 99999.;
1342     }
1343   }
1344
1345   Double_t d0dummy[3]={0.,0.,0.};
1346   AliAODRecoDecayHF aodDplusMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
1347   return aodDplusMC.ImpParXY();
1348
1349 }
1350 //_________________________________________________________________________________________________
1351 Float_t AliAnalysisTaskSEDplus::GetStrangenessWeights(const AliAODRecoDecayHF3Prong* d, TClonesArray* arrayMC, Float_t factor[3]) const {
1352   // Computes weights to adapt strangeness in MC to data
1353
1354   for(Int_t iprong=0;iprong<3;iprong++){
1355     factor[iprong]=1;
1356     AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
1357     Int_t labd= trad->GetLabel();
1358     if(labd>=0){
1359       AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
1360       if(dau){
1361         Int_t labm = dau->GetMother();
1362         if(labm>=0){
1363           AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
1364           if(mot){
1365             if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
1366               if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
1367               else factor[iprong]=1./.6;
1368               //          fNentries->Fill(12);
1369             }
1370             if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
1371               factor[iprong]=1./0.25;
1372               //                  fNentries->Fill(13);
1373             }//if 3122
1374           }//if(mot)
1375         }//if labm>0
1376       }//if(dau)
1377     }//if labd>=0
1378   }//prong loop
1379
1380   Float_t fact=1.;
1381   for(Int_t k=0;k<3;k++)fact=fact*factor[k];
1382   return fact;
1383
1384 }