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