PWG3/vertexingHF/AliAnalysisTaskSEDplus.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //*************************************************************************
  19 // Class AliAnalysisTaskSEDplus
  20 // AliAnalysisTaskSE for the D+ candidates Invariant Mass Histogram and
  21 //comparison of heavy-flavour decay candidates
  22 // to MC truth (kinematics stored in the AOD)
  23 // Authors: Renu Bala, bala@to.infn.it
  24 // F. Prino, prino@to.infn.it
  25 // G. Ortona, ortona@to.infn.it
  26 /////////////////////////////////////////////////////////////
  27
  28 #include <TClonesArray.h>
  29 #include <TNtuple.h>
  30 #include <TCanvas.h>
  31 #include <TList.h>
  32 #include <TString.h>
  33 #include <TH1F.h>
  34 #include <TDatabasePDG.h>
  35
  36 #include "AliAnalysisManager.h"
  37 #include "AliRDHFCutsDplustoKpipi.h"
  38 #include "AliAODHandler.h"
  39 #include "AliAODEvent.h"
  40 #include "AliAODVertex.h"
  41 #include "AliAODTrack.h"
  42 #include "AliAODMCHeader.h"
  43 #include "AliAODMCParticle.h"
  44 #include "AliAODRecoDecayHF3Prong.h"
  45 #include "AliAnalysisVertexingHF.h"
  46 #include "AliAnalysisTaskSE.h"
  47 #include "AliAnalysisTaskSEDplus.h"
  48 #include "AliNormalizationCounter.h"
  49
  50 ClassImp(AliAnalysisTaskSEDplus)
  51
  52
  53 //________________________________________________________________________
  54 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus():
  55 AliAnalysisTaskSE(),
  56   fOutput(0),
  57   fHistNEvents(0),
  58   fPtVsMass(0),
  59   fPtVsMassTC(0),
  60   fYVsPt(0),
  61   fYVsPtTC(0),
  62   fYVsPtSig(0),
  63   fYVsPtSigTC(0),
  64   fNtupleDplus(0),
  65   fUpmasslimit(1.965),
  66   fLowmasslimit(1.765),
  67   fNPtBins(0),
  68   fBinWidth(0.002),
  69   fListCuts(0),
  70   fRDCutsProduction(0),
  71   fRDCutsAnalysis(0),
  72   fCounter(0),
  73   fFillNtuple(kFALSE),
  74   fReadMC(kFALSE),
  75   fUseStrangeness(kFALSE),
  76   fUseBit(kTRUE),
  77   fDoLS(kFALSE)
  78 {
  79    // Default constructor
  80 }
  81
  82 //________________________________________________________________________
  83 AliAnalysisTaskSEDplus::AliAnalysisTaskSEDplus(const char *name,AliRDHFCutsDplustoKpipi *dpluscutsana,AliRDHFCutsDplustoKpipi *dpluscutsprod,Bool_t fillNtuple):
  84 AliAnalysisTaskSE(name),
  85 fOutput(0),
  86 fHistNEvents(0),
  87 fPtVsMass(0),
  88 fPtVsMassTC(0),
  89 fYVsPt(0),
  90 fYVsPtTC(0),
  91 fYVsPtSig(0),
  92 fYVsPtSigTC(0),
  93 fNtupleDplus(0),
  94 fUpmasslimit(1.965),
  95 fLowmasslimit(1.765),
  96 fNPtBins(0),
  97 fBinWidth(0.002),
  98 fListCuts(0),
  99 fRDCutsProduction(dpluscutsprod),
 100 fRDCutsAnalysis(dpluscutsana),
 101 fCounter(0),
 102 fFillNtuple(fillNtuple),
 103 fReadMC(kFALSE),
 104 fUseStrangeness(kFALSE),
 105 fUseBit(kTRUE),
 106 fDoLS(kFALSE)
 107 {
 108   //
 109   // Standrd constructor
 110   //
 111   //Double_t ptlim[5]={0.,2.,3.,5,9999999.};
 112    //SetPtBinLimit(5, ptlim);
 113   SetPtBinLimit(fRDCutsAnalysis->GetNPtBins()+1,fRDCutsAnalysis->GetPtBinLimits());
 114   // Default constructor
 115    // Output slot #1 writes into a TList container
 116   DefineOutput(1,TList::Class());  //My private output
 117  // Output slot #2 writes cut to private output
 118   //  DefineOutput(2,AliRDHFCutsDplustoKpipi::Class());
 119   DefineOutput(2,TList::Class());
 120 // Output slot #3 writes cut to private output
 121   DefineOutput(3,AliNormalizationCounter::Class());
 122
 123   if(fFillNtuple){
 124     // Output slot #4 writes into a TNtuple container
 125     DefineOutput(4,TNtuple::Class());  //My private output
 126   }
 127 }
 128
 129 //________________________________________________________________________
 130 AliAnalysisTaskSEDplus::~AliAnalysisTaskSEDplus()
 131 {
 132   //
 133   // Destructor
 134   //
 135   if (fOutput) {
 136     delete fOutput;
 137     fOutput = 0;
 138   }
 139   if(fHistNEvents){
 140     delete fHistNEvents;
 141     fHistNEvents=0;
 142   }
 143   for(Int_t i=0;i<3;i++){
 144     if(fHistCentrality[i]){delete fHistCentrality[i]; fHistCentrality[i]=0;}
 145   }
 146   for(Int_t i=0;i<3*fNPtBins;i++){
 147     if(fMassHist[i]){ delete fMassHist[i]; fMassHist[i]=0;}
 148     if(fCosPHist[i]){ delete fCosPHist[i]; fCosPHist[i]=0;}
 149     if(fDLenHist[i]){ delete fDLenHist[i]; fDLenHist[i]=0;}
 150     if(fSumd02Hist[i]){ delete fSumd02Hist[i]; fSumd02Hist[i]=0;}
 151     if(fSigVertHist[i]){ delete fSigVertHist[i]; fSigVertHist[i]=0;}
 152     if(fPtMaxHist[i]){ delete fPtMaxHist[i]; fPtMaxHist[i]=0;}
 153     if(fDCAHist[i]){ delete fDCAHist[i]; fDCAHist[i]=0;}
 154     if(fMassHistTC[i]){ delete fMassHistTC[i]; fMassHistTC[i]=0;}
 155     if(fMassHistTCPlus[i]){ delete fMassHistTCPlus[i]; fMassHistTCPlus[i]=0;}
 156     if(fMassHistTCMinus[i]){ delete fMassHistTCMinus[i]; fMassHistTCMinus[i]=0;}
 157
 158     if(fMassHistLS[i]){ delete fMassHistLS[i]; fMassHistLS[i]=0;}
 159     if(fCosPHistLS[i]){ delete fCosPHistLS[i]; fCosPHistLS[i]=0;}
 160     if(fDLenHistLS[i]){ delete fDLenHistLS[i]; fDLenHistLS[i]=0;}
 161     if(fSumd02HistLS[i]){ delete fSumd02HistLS[i]; fSumd02HistLS[i]=0;}
 162     if(fSigVertHistLS[i]){ delete fSigVertHistLS[i]; fSigVertHistLS[i]=0;}
 163     if(fPtMaxHistLS[i]){ delete fPtMaxHistLS[i]; fPtMaxHistLS[i]=0;}
 164     if(fDCAHistLS[i]){ delete fDCAHistLS[i]; fDCAHistLS[i]=0;}
 165     if(fMassHistLSTC[i]){ delete fMassHistLSTC[i]; fMassHistLSTC[i]=0;}
 166   }
 167   if(fPtVsMass){
 168     delete fPtVsMass;
 169     fPtVsMass=0;
 170   }
 171   if(fPtVsMassTC){
 172     delete fPtVsMassTC;
 173     fPtVsMassTC=0;
 174   }
 175   if(fYVsPt){
 176     delete fYVsPt;
 177     fYVsPt=0;
 178   }
 179   if(fYVsPtTC){
 180     delete fYVsPtTC;
 181     fYVsPtTC=0;
 182   }
 183   if(fYVsPtSig){
 184     delete fYVsPtSig;
 185     fYVsPtSig=0;
 186   }
 187   if(fYVsPtSigTC){
 188     delete fYVsPtSigTC;
 189     fYVsPtSigTC=0;
 190   }
 191   if(fNtupleDplus){
 192     delete fNtupleDplus;
 193     fNtupleDplus=0;
 194   }
 195   if (fListCuts) {
 196     delete fListCuts;
 197     fListCuts = 0;
 198   }
 199   if(fRDCutsProduction){
 200     delete fRDCutsProduction;
 201     fRDCutsProduction = 0;
 202   }
 203   if(fRDCutsAnalysis){
 204     delete fRDCutsAnalysis;
 205     fRDCutsAnalysis = 0;
 206   }
 207
 208   if(fCounter){
 209     delete fCounter;
 210     fCounter = 0;
 211   }
 212
 213
 214 }
 215 //_________________________________________________________________
 216 void  AliAnalysisTaskSEDplus::SetMassLimits(Float_t range){
 217   // set invariant mass limits
 218   Float_t bw=GetBinWidth();
 219   fUpmasslimit = 1.865+range;
 220   fLowmasslimit = 1.865-range;
 221   SetBinWidth(bw);
 222 }
 223 //_________________________________________________________________
 224 void  AliAnalysisTaskSEDplus::SetMassLimits(Float_t lowlimit, Float_t uplimit){
 225   // set invariant mass limits
 226   if(uplimit>lowlimit)
 227     {
 228       Float_t bw=GetBinWidth();
 229       fUpmasslimit = lowlimit;
 230       fLowmasslimit = uplimit;
 231       SetBinWidth(bw);
 232     }
 233 }
 234 //________________________________________________________________________
 235 void AliAnalysisTaskSEDplus::SetPtBinLimit(Int_t n, Float_t* lim){
 236   // define pt bins for analysis
 237   if(n>kMaxPtBins){
 238     printf("Max. number of Pt bins = %d\n",kMaxPtBins);
 239     fNPtBins=kMaxPtBins;
 240     fArrayBinLimits[0]=0.;
 241     fArrayBinLimits[1]=2.;
 242     fArrayBinLimits[2]=3.;
 243     fArrayBinLimits[3]=5.;
 244     for(Int_t i=4; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
 245   }else{
 246     fNPtBins=n-1;
 247     fArrayBinLimits[0]=lim[0];
 248     for(Int_t i=1; i<fNPtBins+1; i++)
 249       if(lim[i]>fArrayBinLimits[i-1]){
 250         fArrayBinLimits[i]=lim[i];
 251       }
 252       else {
 253         fArrayBinLimits[i]=fArrayBinLimits[i-1];
 254       }
 255     for(Int_t i=fNPtBins; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
 256   }
 257   if(fDebug > 1){
 258     printf("Number of Pt bins = %d\n",fNPtBins);
 259     for(Int_t i=0; i<fNPtBins; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fArrayBinLimits[i],fArrayBinLimits[i+1]);
 260   }
 261 }
 262 //________________________________________________________________
 263 void AliAnalysisTaskSEDplus::SetBinWidth(Float_t w){
 264   Float_t width=w;
 265   Int_t nbins=(Int_t)((fUpmasslimit-fLowmasslimit)/width+0.5);
 266   Int_t missingbins=4-nbins%4;
 267   nbins=nbins+missingbins;
 268   width=(fUpmasslimit-fLowmasslimit)/nbins;
 269   if(missingbins!=0){
 270     printf("AliAnalysisTaskSEDplus::SetBinWidth: W-bin width of %f will produce histograms not rebinnable by 4. New width set to %f\n",w,width);
 271   }
 272   else{
 273     if(fDebug>1) printf("AliAnalysisTaskSEDplus::SetBinWidth: width set to %f\n",width);
 274   }
 275   fBinWidth=width;
 276 }
 277 //_________________________________________________________________
 278 Double_t  AliAnalysisTaskSEDplus::GetPtBinLimit(Int_t ibin){
 279   // get pt bin limit
 280   if(ibin>fNPtBins)return -1;
 281   return fArrayBinLimits[ibin];
 282 }
 283 //_________________________________________________________________
 284 Int_t AliAnalysisTaskSEDplus::GetNBinsHistos(){
 285   return (Int_t)((fUpmasslimit-fLowmasslimit)/fBinWidth+0.5);
 286 }
 287 //_________________________________________________________________
 288 void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS){
 289   //
 290   //
 291   // Fill the Like Sign histograms
 292   //
 293   printf("started LS\n");
 294   Bool_t fPlotVariables=kTRUE;
 295   //histograms for like sign
 296   Int_t nbins=GetNBinsHistos();;
 297   TH1F *histLSPlus = new TH1F("LSPlus","LSPlus",nbins,fLowmasslimit,fUpmasslimit);
 298   TH1F *histLSMinus = new TH1F("LSMinus","LSMinus",nbins,fLowmasslimit,fUpmasslimit);
 299
 300   Int_t nPosTrks=0,nNegTrks=0;
 301   Int_t nOStriplets = arrayOppositeSign->GetEntriesFast();
 302   Int_t nDplusLS=0;
 303   Int_t nDminusLS=0;
 304   Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
 305   Int_t index=0;
 306
 307   // loop over like sign candidates
 308   for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
 309     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
 310     if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts))continue;
 311     Bool_t unsetvtx=kFALSE;
 312     if(!d->GetOwnPrimaryVtx()) {
 313       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 314       unsetvtx=kTRUE;
 315     }
 316     Double_t ptCand = d->Pt();
 317     Int_t iPtBin = fRDCutsAnalysis->PtBin(ptCand);
 318     if(iPtBin<0)continue;
 319     Int_t sign= d->GetCharge();
 320     if(sign>0){
 321       nPosTrks++;
 322     }else{
 323       nNegTrks++;
 324     }
 325     //    if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)){
 326     fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
 327     Int_t passTightCuts=fRDCutsAnalysis->GetIsSelectedCuts();
 328     if(passTightCuts>0){
 329       Float_t invMass = d->InvMassDplus();
 330       index=GetLSHistoIndex(iPtBin);
 331       fMassHistLS[index+1]->Fill(invMass);
 332       if(sign>0){
 333         histLSPlus->Fill(invMass);
 334         nDplusLS++;
 335       }else{
 336         histLSMinus->Fill(invMass);
 337         nDminusLS++;
 338       }
 339       if(fPlotVariables){
 340         Double_t dlen=d->DecayLength();
 341         Double_t cosp=d->CosPointingAngle();
 342         Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 343         Double_t dca=d->GetDCA();
 344         Double_t sigvert=d->GetSigmaVert();
 345         Double_t ptmax=0;
 346         for(Int_t i=0;i<3;i++){
 347           if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
 348         }
 349         fCosPHistLS[iPtBin]->Fill(cosp);
 350         fDLenHistLS[iPtBin]->Fill(dlen);
 351         fSumd02HistLS[iPtBin]->Fill(sumD02);
 352         fSigVertHistLS[iPtBin]->Fill(sigvert);
 353         fPtMaxHistLS[iPtBin]->Fill(ptmax);
 354         fDCAHistLS[iPtBin]->Fill(dca);
 355       }
 356     }
 357     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 358   }
 359   //wei:normalized to the number of combinations (production)
 360   //wei2:normalized to the number of  LS/OS (production)
 361   //wei3:normalized to the number of  LS/OS (analysis)
 362   //wei4:normalized to the number of combinations (analysis)
 363   Float_t wei2=0;
 364   if(nLikeSign!=0)wei2 = (Float_t)nOStriplets/(Float_t)nLikeSign;
 365   Float_t wei3=0;
 366   if(nDplusLS!=0)wei3 = (Float_t)nDplusOS/(Float_t)nDplusLS;
 367   Float_t weiplus=1.,weiminus=1.;
 368   Float_t wei4plus=1.,wei4minus=1.;
 369   //wei* should be automatically protected, since to get a triplet there must be at least 3 good tracks in the event
 370   if(nPosTrks>2)weiplus=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.);
 371   if(nDplusLS>2)wei4plus=3.*(Float_t)nDminusLS/((Float_t)nDplusLS-2.);
 372   if(nNegTrks>2)weiminus=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.);
 373   if(nDminusLS>2)wei4minus=3.*(Float_t)nDplusLS/((Float_t)nDminusLS-2.);
 374
 375   fMassHistLS[index]->Add(histLSPlus,weiplus);
 376   fMassHistLS[index]->Add(histLSMinus,weiminus);
 377   fMassHistLS[index+2]->Add(histLSPlus,wei2);
 378   fMassHistLS[index+2]->Add(histLSMinus,wei2);
 379   fMassHistLS[index+3]->Add(histLSPlus,wei3);
 380   fMassHistLS[index+3]->Add(histLSMinus,wei3);
 381   fMassHistLS[index+4]->Add(histLSPlus,wei4plus);
 382   fMassHistLS[index+4]->Add(histLSMinus,wei4minus);
 383
 384   delete histLSPlus;histLSPlus=0;
 385   delete histLSMinus;histLSMinus=0;
 386
 387   if(fDebug>1) printf("LS analysis terminated\n");
 388 }
 389
 390 //__________________________________________
 391 void AliAnalysisTaskSEDplus::Init(){
 392   //
 393   // Initialization
 394   //
 395   if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
 396
 397   //PostData(2,fRDCutsloose);//we should then put those cuts in a tlist if we have more than 1
 398   fListCuts=new TList();
 399   AliRDHFCutsDplustoKpipi *production = new AliRDHFCutsDplustoKpipi(*fRDCutsProduction);
 400   production->SetName("ProductionCuts");
 401   AliRDHFCutsDplustoKpipi *analysis = new AliRDHFCutsDplustoKpipi(*fRDCutsAnalysis);
 402   analysis->SetName("AnalysisCuts");
 403
 404   fListCuts->Add(production);
 405   fListCuts->Add(analysis);
 406   PostData(2,fListCuts);
 407
 408   return;
 409 }
 410
 411 //________________________________________________________________________
 412 void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
 413 {
 414   // Create the output container
 415   //
 416   if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");
 417
 418   // Several histograms are more conveniently managed in a TList
 419   fOutput = new TList();
 420   fOutput->SetOwner();
 421   fOutput->SetName("OutputHistos");
 422
 423   TString hisname;
 424   Int_t index=0;
 425   Int_t indexLS=0;
 426   Int_t nbins=GetNBinsHistos();
 427   fHistCentrality[0]=new TH1F("centrality","centrality",100,0.5,30000.5);
 428   fHistCentrality[1]=new TH1F("centrality(selectedCent)","centrality(selectedCent)",100,0.5,30000.5);
 429   fHistCentrality[2]=new TH1F("centrality(OutofCent)","centrality(OutofCent)",100,0.5,30000.5);
 430   for(Int_t i=0;i<3;i++){
 431     fHistCentrality[i]->Sumw2();
 432     fOutput->Add(fHistCentrality[i]);
 433   }
 434   for(Int_t i=0;i<fNPtBins;i++){
 435
 436     index=GetHistoIndex(i);
 437     indexLS=GetLSHistoIndex(i);
 438
 439     hisname.Form("hMassPt%d",i);
 440     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 441     fMassHist[index]->Sumw2();
 442     hisname.Form("hCosPAllPt%d",i);
 443     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
 444     fCosPHist[index]->Sumw2();
 445     hisname.Form("hDLenAllPt%d",i);
 446     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
 447     fDLenHist[index]->Sumw2();
 448     hisname.Form("hSumd02AllPt%d",i);
 449     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
 450     fSumd02Hist[index]->Sumw2();
 451     hisname.Form("hSigVertAllPt%d",i);
 452     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 453     fSigVertHist[index]->Sumw2();
 454     hisname.Form("hPtMaxAllPt%d",i);
 455     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 456     fPtMaxHist[index]->Sumw2();
 457
 458     hisname.Form("hDCAAllPt%d",i);
 459     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 460     fDCAHist[index]->Sumw2();
 461
 462
 463
 464     hisname.Form("hMassPt%dTC",i);
 465     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 466     fMassHistTC[index]->Sumw2();
 467     hisname.Form("hMassPt%dTCPlus",i);
 468     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 469     fMassHistTCPlus[index]->Sumw2();
 470     hisname.Form("hMassPt%dTCMinus",i);
 471     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 472     fMassHistTCMinus[index]->Sumw2();
 473
 474
 475
 476
 477     hisname.Form("hCosPAllPt%dLS",i);
 478     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
 479     fCosPHistLS[index]->Sumw2();
 480     hisname.Form("hDLenAllPt%dLS",i);
 481     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
 482     fDLenHistLS[index]->Sumw2();
 483     hisname.Form("hSumd02AllPt%dLS",i);
 484     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
 485     fSumd02HistLS[index]->Sumw2();
 486     hisname.Form("hSigVertAllPt%dLS",i);
 487     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 488     fSigVertHistLS[index]->Sumw2();
 489     hisname.Form("hPtMaxAllPt%dLS",i);
 490     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 491     fPtMaxHistLS[index]->Sumw2();
 492
 493     hisname.Form("hDCAAllPt%dLS",i);
 494     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 495     fDCAHistLS[index]->Sumw2();
 496
 497     hisname.Form("hLSPt%dLC",i);
 498     fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 499     fMassHistLS[indexLS]->Sumw2();
 500
 501     hisname.Form("hLSPt%dTC",i);
 502     fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 503     fMassHistLSTC[indexLS]->Sumw2();
 504
 505
 506
 507     index=GetSignalHistoIndex(i);
 508     indexLS++;
 509     hisname.Form("hSigPt%d",i);
 510     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 511     fMassHist[index]->Sumw2();
 512     hisname.Form("hCosPSigPt%d",i);
 513     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
 514     fCosPHist[index]->Sumw2();
 515     hisname.Form("hDLenSigPt%d",i);
 516     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
 517     fDLenHist[index]->Sumw2();
 518     hisname.Form("hSumd02SigPt%d",i);
 519     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
 520     fSumd02Hist[index]->Sumw2();
 521     hisname.Form("hSigVertSigPt%d",i);
 522     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 523     fSigVertHist[index]->Sumw2();
 524     hisname.Form("hPtMaxSigPt%d",i);
 525     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 526     fPtMaxHist[index]->Sumw2();
 527
 528     hisname.Form("hDCASigPt%d",i);
 529     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 530     fDCAHist[index]->Sumw2();
 531
 532
 533     hisname.Form("hSigPt%dTC",i);
 534     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 535     fMassHistTC[index]->Sumw2();
 536     hisname.Form("hSigPt%dTCPlus",i);
 537     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 538     fMassHistTCPlus[index]->Sumw2();
 539     hisname.Form("hSigPt%dTCMinus",i);
 540     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 541     fMassHistTCMinus[index]->Sumw2();
 542
 543     hisname.Form("hLSPt%dLCnw",i);
 544     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 545     fMassHistLS[indexLS]->Sumw2();
 546     hisname.Form("hLSPt%dTCnw",i);
 547     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 548     fMassHistLSTC[indexLS]->Sumw2();
 549
 550
 551
 552     hisname.Form("hCosPSigPt%dLS",i);
 553     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
 554     fCosPHistLS[index]->Sumw2();
 555     hisname.Form("hDLenSigPt%dLS",i);
 556     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
 557     fDLenHistLS[index]->Sumw2();
 558     hisname.Form("hSumd02SigPt%dLS",i);
 559     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
 560     fSumd02HistLS[index]->Sumw2();
 561     hisname.Form("hSigVertSigPt%dLS",i);
 562     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 563     fSigVertHistLS[index]->Sumw2();
 564     hisname.Form("hPtMaxSigPt%dLS",i);
 565     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 566     fPtMaxHistLS[index]->Sumw2();
 567
 568     hisname.Form("hDCASigPt%dLS",i);
 569     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 570     fDCAHistLS[index]->Sumw2();
 571
 572
 573
 574     index=GetBackgroundHistoIndex(i);
 575     indexLS++;
 576     hisname.Form("hBkgPt%d",i);
 577     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 578     fMassHist[index]->Sumw2();
 579     hisname.Form("hCosPBkgPt%d",i);
 580     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
 581     fCosPHist[index]->Sumw2();
 582     hisname.Form("hDLenBkgPt%d",i);
 583     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
 584     fDLenHist[index]->Sumw2();
 585     hisname.Form("hSumd02BkgPt%d",i);
 586     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
 587     fSumd02Hist[index]->Sumw2();
 588     hisname.Form("hSigVertBkgPt%d",i);
 589     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 590     fSigVertHist[index]->Sumw2();
 591     hisname.Form("hPtMaxBkgPt%d",i);
 592     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 593     fPtMaxHist[index]->Sumw2();
 594
 595     hisname.Form("hDCABkgPt%d",i);
 596     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 597     fDCAHist[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     hisname.Form("hLSPt%dLCntrip",i);
 611     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 612     fMassHistLS[indexLS]->Sumw2();
 613     hisname.Form("hLSPt%dTCntrip",i);
 614     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 615     fMassHistLSTC[indexLS]->Sumw2();
 616
 617
 618     hisname.Form("hCosPBkgPt%dLS",i);
 619     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
 620     fCosPHistLS[index]->Sumw2();
 621     hisname.Form("hDLenBkgPt%dLS",i);
 622     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
 623     fDLenHistLS[index]->Sumw2();
 624     hisname.Form("hSumd02BkgPt%dLS",i);
 625     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
 626     fSumd02HistLS[index]->Sumw2();
 627     hisname.Form("hSigVertBkgPt%dLS",i);
 628     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 629     fSigVertHistLS[index]->Sumw2();
 630     hisname.Form("hPtMaxBkgPt%dLS",i);
 631     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 632     fPtMaxHistLS[index]->Sumw2();
 633     hisname.Form("hDCABkgPt%dLS",i);
 634     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 635     fDCAHistLS[index]->Sumw2();
 636
 637
 638     indexLS++;
 639     hisname.Form("hLSPt%dLCntripsinglecut",i);
 640     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 641     fMassHistLS[indexLS]->Sumw2();
 642     hisname.Form("hLSPt%dTCntripsinglecut",i);
 643     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 644     fMassHistLSTC[indexLS]->Sumw2();
 645
 646     indexLS++;
 647     hisname.Form("hLSPt%dLCspc",i);
 648     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 649     fMassHistLS[indexLS]->Sumw2();
 650     hisname.Form("hLSPt%dTCspc",i);
 651     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 652     fMassHistLSTC[indexLS]->Sumw2();
 653   }
 654
 655   for(Int_t i=0; i<3*fNPtBins; i++){
 656     fOutput->Add(fMassHist[i]);
 657      fOutput->Add(fCosPHist[i]);
 658     fOutput->Add(fDLenHist[i]);
 659     fOutput->Add(fSumd02Hist[i]);
 660     fOutput->Add(fSigVertHist[i]);
 661     fOutput->Add(fPtMaxHist[i]);
 662     fOutput->Add(fDCAHist[i]);
 663     fOutput->Add(fMassHistTC[i]);
 664     fOutput->Add(fMassHistTCPlus[i]);
 665     fOutput->Add(fMassHistTCMinus[i]);
 666   }
 667   for(Int_t i=0; i<3*fNPtBins&&fDoLS; i++){
 668     fOutput->Add(fCosPHistLS[i]);
 669     fOutput->Add(fDLenHistLS[i]);
 670     fOutput->Add(fSumd02HistLS[i]);
 671     fOutput->Add(fSigVertHistLS[i]);
 672     fOutput->Add(fPtMaxHistLS[i]);
 673     fOutput->Add(fDCAHistLS[i]);
 674   }
 675   for(Int_t i=0; i<5*fNPtBins&&fDoLS; i++){
 676     fOutput->Add(fMassHistLS[i]);
 677     fOutput->Add(fMassHistLSTC[i]);
 678   }
 679
 680
 681   fHistNEvents = new TH1F("fHistNEvents", "number of events ",9,-0.5,8.5);
 682   fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 683   fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with good vertex");
 684   fHistNEvents->GetXaxis()->SetBinLabel(3,"nEvents with PbPb HM trigger");
 685   fHistNEvents->GetXaxis()->SetBinLabel(4,"no. of Rejected pileup events");
 686   fHistNEvents->GetXaxis()->SetBinLabel(5,"no. of candidate");
 687   fHistNEvents->GetXaxis()->SetBinLabel(6,"no. of D+ after loose cuts");
 688   fHistNEvents->GetXaxis()->SetBinLabel(7,"no. of D+ after tight cuts");
 689   fHistNEvents->GetXaxis()->SetBinLabel(8,"no. of out centrality events");
 690   fHistNEvents->GetXaxis()->SetBinLabel(9,"no. of cand wo bitmask");
 691
 692   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
 693
 694   fHistNEvents->Sumw2();
 695   fHistNEvents->SetMinimum(0);
 696   fOutput->Add(fHistNEvents);
 697
 698   fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
 699   fPtVsMassTC=new TH2F("hPtVsMassTC","PtVsMass (analysis cuts)",nbins,fLowmasslimit,fUpmasslimit,200,0.,20.);
 700   fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
 701   fYVsPtTC=new TH2F("hYVsPtTC","YvsPt (analysis cuts)",40,0.,20.,80,-2.,2.);
 702   fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
 703   fYVsPtSigTC=new TH2F("hYVsPtSigTC","YvsPt (MC, only Sig, analysis cuts)",40,0.,20.,80,-2.,2.);
 704
 705   fOutput->Add(fPtVsMass);
 706   fOutput->Add(fPtVsMassTC);
 707   fOutput->Add(fYVsPt);
 708   fOutput->Add(fYVsPtTC);
 709   fOutput->Add(fYVsPtSig);
 710   fOutput->Add(fYVsPtSigTC);
 711
 712
 713   //Counter for Normalization
 714   TString normName="NormalizationCounter";
 715   AliAnalysisDataContainer *cont = GetOutputSlot(3)->GetContainer();
 716   if(cont)normName=(TString)cont->GetName();
 717   fCounter = new AliNormalizationCounter(normName.Data());
 718   fCounter->SetRejectPileUp(fRDCutsProduction->GetOptPileUp());
 719
 720   if(fFillNtuple){
 721     OpenFile(4); // 4 is the slot number of the ntuple
 722
 723     fNtupleDplus = new TNtuple("fNtupleDplus","D +","pdg:Px:Py:Pz:PtTrue:VxTrue:VyTrue:VzTrue:Ptpi:PtK:Ptpi2:PtRec:PointingAngle:DecLeng:VxRec:VyRec:VzRec:InvMass:sigvert:d0Pi:d0K:d0Pi2:dca:d0square");
 724
 725   }
 726
 727   return;
 728 }
 729
 730 //________________________________________________________________________
 731 void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
 732 {
 733   // Execute analysis for current event:
 734   // heavy flavor candidates association to MC truth
 735
 736   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 737
 738
 739
 740   TClonesArray *array3Prong = 0;
 741   TClonesArray *arrayLikeSign =0;
 742   if(!aod && AODEvent() && IsStandardAOD()) {
 743     // In case there is an AOD handler writing a standard AOD, use the AOD
 744     // event in memory rather than the input (ESD) event.
 745     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 746      // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 747      // have to taken from the AOD event hold by the AliAODExtension
 748     AliAODHandler* aodHandler = (AliAODHandler*)
 749       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 750     if(aodHandler->GetExtensions()) {
 751       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 752       AliAODEvent *aodFromExt = ext->GetAOD();
 753       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 754       arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
 755     }
 756   } else if(aod) {
 757     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 758     arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
 759   }
 760
 761   if(!aod || !array3Prong) {
 762     printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
 763     return;
 764   }
 765   if(!arrayLikeSign && fDoLS) {
 766     printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
 767     return;
 768   }
 769
 770
 771   // fix for temporary bug in ESDfilter
 772   // the AODs with null vertex pointer didn't pass the PhysSel
 773   if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
 774   fCounter->StoreEvent(aod,fReadMC);
 775   fHistNEvents->Fill(0); // count event
 776
 777   Bool_t isEvSel=fRDCutsAnalysis->IsEventSelected(aod);
 778   Float_t centrality=aod->GetNTracks();//fRDCutsAnalysis->GetCentrality(aod);
 779   fHistCentrality[0]->Fill(centrality);
 780   // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
 781   TString trigclass=aod->GetFiredTriggerClasses();
 782   if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")||trigclass.Contains("C0SMH-B-NOPF-ALL")) fHistNEvents->Fill(2);
 783   if(fRDCutsAnalysis->GetWhyRejection()==1)fHistNEvents->Fill(3);
 784   if(fRDCutsAnalysis->GetWhyRejection()==2){fHistNEvents->Fill(7);fHistCentrality[2]->Fill(centrality);}
 785
 786   // Post the data already here
 787   PostData(1,fOutput);
 788   if(!isEvSel)return;
 789
 790   fHistCentrality[1]->Fill(centrality);
 791   fHistNEvents->Fill(1);
 792
 793   TClonesArray *arrayMC=0;
 794   AliAODMCHeader *mcHeader=0;
 795
 796   // AOD primary vertex
 797   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 798   //    vtx1->Print();
 799    TString primTitle = vtx1->GetTitle();
 800    //if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0)fHistNEvents->Fill(2);
 801
 802   // load MC particles
 803   if(fReadMC){
 804
 805     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 806     if(!arrayMC) {
 807       printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
 808       return;
 809     }
 810
 811   // load MC header
 812     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 813     if(!mcHeader) {
 814     printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
 815     return;
 816     }
 817   }
 818
 819   Int_t n3Prong = array3Prong->GetEntriesFast();
 820   printf("Number of D+->Kpipi: %d and of tracks: %d\n",n3Prong,aod->GetNumberOfTracks());
 821
 822
 823   Int_t nOS=0;
 824   Int_t index;
 825   Int_t pdgDgDplustoKpipi[3]={321,211,211};
 826   // 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
 827   //Double_t *cutsDplus = new (Double_t*)fRDCuts->GetCuts();
 828   Int_t nSelectedloose=0,nSelectedtight=0;
 829   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 830     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 831     fHistNEvents->Fill(4);
 832     if(fUseBit && !d->HasSelectionBit(AliRDHFCuts::kDplusCuts)){
 833       fHistNEvents->Fill(8);
 834       continue;
 835     }
 836     Bool_t unsetvtx=kFALSE;
 837     if(!d->GetOwnPrimaryVtx()){
 838       d->SetOwnPrimaryVtx(vtx1);
 839       unsetvtx=kTRUE;
 840     }
 841
 842     if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate,aod)) {
 843
 844
 845
 846
 847       Int_t iPtBin = -1;
 848       Double_t ptCand = d->Pt();
 849
 850       for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
 851         if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
 852       }
 853
 854       Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate,aod);
 855
 856
 857       Int_t labDp=-1;
 858       Float_t deltaPx=0.;
 859       Float_t deltaPy=0.;
 860       Float_t deltaPz=0.;
 861       Float_t truePt=0.;
 862       Float_t xDecay=0.;
 863       Float_t yDecay=0.;
 864       Float_t zDecay=0.;
 865       Float_t pdgCode=-2;
 866       if(fReadMC){
 867         labDp = d->MatchToMC(411,arrayMC,3,pdgDgDplustoKpipi);
 868         if(labDp>=0){
 869           AliAODMCParticle *partDp = (AliAODMCParticle*)arrayMC->At(labDp);
 870           AliAODMCParticle *dg0 = (AliAODMCParticle*)arrayMC->At(partDp->GetDaughter(0));
 871           deltaPx=partDp->Px()-d->Px();
 872           deltaPy=partDp->Py()-d->Py();
 873           deltaPz=partDp->Pz()-d->Pz();
 874           truePt=partDp->Pt();
 875           xDecay=dg0->Xv();
 876           yDecay=dg0->Yv();
 877           zDecay=dg0->Zv();
 878           pdgCode=TMath::Abs(partDp->GetPdgCode());
 879         }else{
 880           pdgCode=-1;
 881         }
 882       }
 883       Double_t invMass=d->InvMassDplus();
 884       Double_t rapid=d->YDplus();
 885       fYVsPt->Fill(ptCand,rapid);
 886       if(passTightCuts) fYVsPtTC->Fill(ptCand,rapid);
 887       Bool_t isFidAcc=fRDCutsAnalysis->IsInFiducialAcceptance(ptCand,rapid);
 888       if(isFidAcc){
 889         fPtVsMass->Fill(invMass,ptCand);
 890         if(passTightCuts) fPtVsMassTC->Fill(invMass,ptCand);
 891       }
 892       Float_t tmp[24];
 893       if(fFillNtuple){
 894         tmp[0]=pdgCode;
 895         tmp[1]=deltaPx;
 896         tmp[2]=deltaPy;
 897         tmp[3]=deltaPz;
 898         tmp[4]=truePt;
 899         tmp[5]=xDecay;
 900         tmp[6]=yDecay;
 901         tmp[7]=zDecay;
 902         tmp[8]=d->PtProng(0);
 903         tmp[9]=d->PtProng(1);
 904         tmp[10]=d->PtProng(2);
 905         tmp[11]=d->Pt();
 906         tmp[12]=d->CosPointingAngle();
 907         tmp[13]=d->DecayLength();
 908         tmp[14]=d->Xv();
 909         tmp[15]=d->Yv();
 910         tmp[16]=d->Zv();
 911         tmp[17]=d->InvMassDplus();
 912         tmp[18]=d->GetSigmaVert();
 913         tmp[19]=d->Getd0Prong(0);
 914         tmp[20]=d->Getd0Prong(1);
 915         tmp[21]=d->Getd0Prong(2);
 916         tmp[22]=d->GetDCA();
 917         tmp[23]=d->Prodd0d0();
 918         fNtupleDplus->Fill(tmp);
 919         PostData(4,fNtupleDplus);
 920       }
 921       Double_t dlen=d->DecayLength();
 922       Double_t cosp=d->CosPointingAngle();
 923       Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 924       Double_t dca=d->GetDCA();
 925       Double_t sigvert=d->GetSigmaVert();
 926       Double_t ptmax=0;
 927       for(Int_t i=0;i<3;i++){
 928         if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
 929       }
 930       if(iPtBin>=0){
 931
 932         index=GetHistoIndex(iPtBin);
 933         if(isFidAcc){
 934           fHistNEvents->Fill(5);
 935           nSelectedloose++;
 936           fMassHist[index]->Fill(invMass);
 937           fCosPHist[index]->Fill(cosp);
 938           fDLenHist[index]->Fill(dlen);
 939           fSumd02Hist[index]->Fill(sumD02);
 940           fSigVertHist[index]->Fill(sigvert);
 941           fPtMaxHist[index]->Fill(ptmax);
 942           fDCAHist[index]->Fill(dca);
 943
 944           if(passTightCuts){ fHistNEvents->Fill(6);
 945             nSelectedtight++;
 946             fMassHistTC[index]->Fill(invMass);
 947             if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
 948             else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
 949           }
 950         }
 951
 952         if(fReadMC){
 953           if(labDp>=0) {
 954             index=GetSignalHistoIndex(iPtBin);
 955             if(isFidAcc){
 956               Float_t factor[3]={1.,1.,1.};
 957               if(fUseStrangeness){
 958                 for(Int_t iprong=0;iprong<3;iprong++){
 959                   AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
 960                   Int_t labd= trad->GetLabel();
 961                   if(labd>=0){
 962                     AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
 963                     if(dau){
 964                       Int_t labm = dau->GetMother();
 965                       if(labm>=0){
 966                         AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
 967                         if(mot){
 968                           if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
 969                             if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
 970                             else factor[iprong]=1./.6;
 971                             //          fNentries->Fill(12);
 972                           }
 973                           if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
 974                             factor[iprong]=1./0.25;
 975                             //            fNentries->Fill(13);
 976                           }//if 3122
 977                         }//if(mot)
 978                       }//if labm>0
 979                     }//if(dau)
 980                   }//if labd>=0
 981                 }//prong loop
 982               }
 983               Float_t fact=1.;for(Int_t k=0;k<3;k++)fact=fact*factor[k];
 984               fMassHist[index]->Fill(invMass);
 985               fCosPHist[index]->Fill(cosp,fact);
 986               fDLenHist[index]->Fill(dlen,fact);
 987               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];
 988               fSumd02Hist[index]->Fill(sumd02s);
 989               fSigVertHist[index]->Fill(sigvert,fact);
 990               fPtMaxHist[index]->Fill(ptmax,fact);
 991               fDCAHist[index]->Fill(dca,fact);
 992               if(passTightCuts){
 993                 fMassHistTC[index]->Fill(invMass);
 994                 if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
 995                 else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
 996               }
 997             }
 998             fYVsPtSig->Fill(ptCand,rapid);
 999             if(passTightCuts) fYVsPtSigTC->Fill(ptCand,rapid);
1000           }else{
1001             index=GetBackgroundHistoIndex(iPtBin);
1002             if(isFidAcc){
1003               Float_t factor[3]={1.,1.,1.};
1004               if(fUseStrangeness){
1005                 for(Int_t iprong=0;iprong<3;iprong++){
1006                   AliAODTrack *trad = (AliAODTrack*)d->GetDaughter(iprong);
1007                   Int_t labd= trad->GetLabel();
1008                   if(labd>=0){
1009                     AliAODMCParticle *dau = (AliAODMCParticle*)arrayMC->At(labd);
1010                     if(dau){
1011                       Int_t labm = dau->GetMother();
1012                       if(labm>=0){
1013                         AliAODMCParticle *mot = (AliAODMCParticle*)arrayMC->At(labm);
1014                         if(mot){
1015                           if(TMath::Abs(mot->GetPdgCode())==310 || TMath::Abs(mot->GetPdgCode())==130 || TMath::Abs(mot->GetPdgCode())==321){ //K0_S, K0_L, K^+-
1016                             if(d->PtProng(iprong)<=1)factor[iprong]=1./.7;
1017                             else factor[iprong]=1./.6;
1018                             //          fNentries->Fill(12);
1019                           }
1020                           if(TMath::Abs(mot->GetPdgCode())==3122) { //Lambda
1021                             factor[iprong]=1./0.25;
1022                             //            fNentries->Fill(13);
1023                           }//if 3122
1024                         }//if(mot)
1025                       }//if labm>0
1026                     }//if(dau)
1027                   }//if labd>=0
1028                 }//prong loop
1029               }
1030               Float_t fact=1.;for(Int_t k=0;k<3;k++)fact=fact*factor[k];
1031               fMassHist[index]->Fill(invMass);
1032               fCosPHist[index]->Fill(cosp,fact);
1033               fDLenHist[index]->Fill(dlen,fact);
1034               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];
1035               fSumd02Hist[index]->Fill(sumd02s);
1036               fSigVertHist[index]->Fill(sigvert,fact);
1037               fPtMaxHist[index]->Fill(ptmax,fact);
1038               fDCAHist[index]->Fill(dca,fact);
1039               if(passTightCuts){
1040                 fMassHistTC[index]->Fill(invMass);
1041                 if(d->GetCharge()>0) fMassHistTCPlus[index]->Fill(invMass);
1042                 else if(d->GetCharge()<0) fMassHistTCMinus[index]->Fill(invMass);
1043               }
1044             }
1045           }
1046
1047           }
1048
1049     }
1050
1051     }
1052     if(unsetvtx) d->UnsetOwnPrimaryVtx();
1053   }
1054   fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
1055   fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
1056
1057   //start LS analysis
1058   if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
1059
1060   PostData(1,fOutput);
1061   PostData(2,fListCuts);
1062   PostData(3,fCounter);
1063   return;
1064 }
1065
1066
1067
1068 //________________________________________________________________________
1069 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
1070 {
1071   // Terminate analysis
1072   //
1073   if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
1074
1075   fOutput = dynamic_cast<TList*> (GetOutputData(1));
1076   if (!fOutput) {
1077     printf("ERROR: fOutput not available\n");
1078     return;
1079   }
1080   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
1081
1082   TString hisname;
1083   Int_t index=0;
1084
1085   for(Int_t i=0;i<fNPtBins;i++){
1086     index=GetHistoIndex(i);
1087
1088     hisname.Form("hMassPt%dTC",i);
1089     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1090   }
1091
1092   TCanvas *c1=new TCanvas("c1","D+ invariant mass distribution",500,500);
1093   c1->cd();
1094   TH1F *hMassPt=(TH1F*)fOutput->FindObject("hMassPt3TC");
1095   hMassPt->SetLineColor(kBlue);
1096   hMassPt->SetXTitle("M[GeV/c^{2}]");
1097   hMassPt->Draw();
1098
1099   return;
1100 }