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*fNPtBins;i++){
 140     if(fMassHist[i]){ delete fMassHist[i]; fMassHist[i]=0;}
 141     if(fCosPHist[i]){ delete fCosPHist[i]; fCosPHist[i]=0;}
 142     if(fDLenHist[i]){ delete fDLenHist[i]; fDLenHist[i]=0;}
 143     if(fSumd02Hist[i]){ delete fSumd02Hist[i]; fSumd02Hist[i]=0;}
 144     if(fSigVertHist[i]){ delete fSigVertHist[i]; fSigVertHist[i]=0;}
 145     if(fPtMaxHist[i]){ delete fPtMaxHist[i]; fPtMaxHist[i]=0;}
 146     if(fDCAHist[i]){ delete fDCAHist[i]; fDCAHist[i]=0;}
 147     if(fMassHistTC[i]){ delete fMassHistTC[i]; fMassHistTC[i]=0;}
 148     if(fMassHistTCPlus[i]){ delete fMassHistTCPlus[i]; fMassHistTCPlus[i]=0;}
 149     if(fMassHistTCMinus[i]){ delete fMassHistTCMinus[i]; fMassHistTCMinus[i]=0;}
 150
 151     if(fMassHistLS[i]){ delete fMassHistLS[i]; fMassHistLS[i]=0;}
 152     if(fCosPHistLS[i]){ delete fCosPHistLS[i]; fCosPHistLS[i]=0;}
 153     if(fDLenHistLS[i]){ delete fDLenHistLS[i]; fDLenHistLS[i]=0;}
 154     if(fSumd02HistLS[i]){ delete fSumd02HistLS[i]; fSumd02HistLS[i]=0;}
 155     if(fSigVertHistLS[i]){ delete fSigVertHistLS[i]; fSigVertHistLS[i]=0;}
 156     if(fPtMaxHistLS[i]){ delete fPtMaxHistLS[i]; fPtMaxHistLS[i]=0;}
 157     if(fDCAHistLS[i]){ delete fDCAHistLS[i]; fDCAHistLS[i]=0;}
 158     if(fMassHistLSTC[i]){ delete fMassHistLSTC[i]; fMassHistLSTC[i]=0;}
 159   }
 160   if(fPtVsMass){
 161     delete fPtVsMass;
 162     fPtVsMass=0;
 163   }
 164   if(fPtVsMassTC){
 165     delete fPtVsMassTC;
 166     fPtVsMassTC=0;
 167   }
 168   if(fYVsPt){
 169     delete fYVsPt;
 170     fYVsPt=0;
 171   }
 172   if(fYVsPtTC){
 173     delete fYVsPtTC;
 174     fYVsPtTC=0;
 175   }
 176   if(fYVsPtSig){
 177     delete fYVsPtSig;
 178     fYVsPtSig=0;
 179   }
 180   if(fYVsPtSigTC){
 181     delete fYVsPtSigTC;
 182     fYVsPtSigTC=0;
 183   }
 184   if(fNtupleDplus){
 185     delete fNtupleDplus;
 186     fNtupleDplus=0;
 187   }
 188   if (fListCuts) {
 189     delete fListCuts;
 190     fListCuts = 0;
 191   }
 192   if(fRDCutsProduction){
 193     delete fRDCutsProduction;
 194     fRDCutsProduction = 0;
 195   }
 196   if(fRDCutsAnalysis){
 197     delete fRDCutsAnalysis;
 198     fRDCutsAnalysis = 0;
 199   }
 200
 201   if(fCounter){
 202     delete fCounter;
 203     fCounter = 0;
 204   }
 205
 206
 207 }
 208 //_________________________________________________________________
 209 void  AliAnalysisTaskSEDplus::SetMassLimits(Float_t range){
 210   // set invariant mass limits
 211   Float_t bw=GetBinWidth();
 212   fUpmasslimit = 1.865+range;
 213   fLowmasslimit = 1.865-range;
 214   SetBinWidth(bw);
 215 }
 216 //_________________________________________________________________
 217 void  AliAnalysisTaskSEDplus::SetMassLimits(Float_t lowlimit, Float_t uplimit){
 218   // set invariant mass limits
 219   if(uplimit>lowlimit)
 220     {
 221       Float_t bw=GetBinWidth();
 222       fUpmasslimit = lowlimit;
 223       fLowmasslimit = uplimit;
 224       SetBinWidth(bw);
 225     }
 226 }
 227 //________________________________________________________________________
 228 void AliAnalysisTaskSEDplus::SetPtBinLimit(Int_t n, Float_t* lim){
 229   // define pt bins for analysis
 230   if(n>kMaxPtBins){
 231     printf("Max. number of Pt bins = %d\n",kMaxPtBins);
 232     fNPtBins=kMaxPtBins;
 233     fArrayBinLimits[0]=0.;
 234     fArrayBinLimits[1]=2.;
 235     fArrayBinLimits[2]=3.;
 236     fArrayBinLimits[3]=5.;
 237     for(Int_t i=4; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
 238   }else{
 239     fNPtBins=n-1;
 240     fArrayBinLimits[0]=lim[0];
 241     for(Int_t i=1; i<fNPtBins+1; i++)
 242       if(lim[i]>fArrayBinLimits[i-1]){
 243         fArrayBinLimits[i]=lim[i];
 244       }
 245       else {
 246         fArrayBinLimits[i]=fArrayBinLimits[i-1];
 247       }
 248     for(Int_t i=fNPtBins; i<kMaxPtBins+1; i++) fArrayBinLimits[i]=99999999.;
 249   }
 250   if(fDebug > 1){
 251     printf("Number of Pt bins = %d\n",fNPtBins);
 252     for(Int_t i=0; i<fNPtBins+1; i++) printf(" Bin%d = %8.2f-%8.2f\n",i,fArrayBinLimits[i],fArrayBinLimits[i+1]);
 253   }
 254 }
 255 //________________________________________________________________
 256 void AliAnalysisTaskSEDplus::SetBinWidth(Float_t w){
 257   Float_t width=w;
 258   Int_t nbins=(Int_t)((fUpmasslimit-fLowmasslimit)/width+0.5);
 259   Int_t missingbins=4-nbins%4;
 260   nbins=nbins+missingbins;
 261   width=(fUpmasslimit-fLowmasslimit)/nbins;
 262   if(missingbins!=0){
 263     printf("AliAnalysisTaskSEDplus::SetBinWidth: W-bin width of %f will produce histograms not rebinnable by 4. New width set to %f\n",w,width);
 264   }
 265   else{
 266     if(fDebug>1) printf("AliAnalysisTaskSEDplus::SetBinWidth: width set to %f\n",width);
 267   }
 268   fBinWidth=width;
 269 }
 270 //_________________________________________________________________
 271 Double_t  AliAnalysisTaskSEDplus::GetPtBinLimit(Int_t ibin){
 272   // get pt bin limit
 273   if(ibin>fNPtBins)return -1;
 274   return fArrayBinLimits[ibin];
 275 }
 276 //_________________________________________________________________
 277 Int_t AliAnalysisTaskSEDplus::GetNBinsHistos(){
 278   return (Int_t)((fUpmasslimit-fLowmasslimit)/fBinWidth+0.5);
 279 }
 280 //_________________________________________________________________
 281 void AliAnalysisTaskSEDplus::LSAnalysis(TClonesArray *arrayOppositeSign,TClonesArray *arrayLikeSign,AliAODEvent *aod,AliAODVertex *vtx1, Int_t nDplusOS){
 282   //
 283   //
 284   // Fill the Like Sign histograms
 285   //
 286
 287   //count pos/neg tracks
 288   Int_t nPosTrks=0,nNegTrks=0;
 289  //counter for particles passing single particle cuts
 290   Int_t nspcplus=0;
 291   Int_t nspcminus=0;
 292
 293   for(Int_t it=0;it<aod->GetNumberOfTracks();it++) {
 294     AliAODTrack *track = aod->GetTrack(it);
 295     if(track->Charge()>0){
 296       nPosTrks++;
 297       if(track->Pt()>=0.4){
 298         nspcplus++;
 299       }
 300     }
 301     if(track->Charge()<0)
 302       {
 303         nNegTrks++;
 304         if(track->Pt()>=0.4){
 305           nspcminus++;
 306         }
 307       }
 308   }
 309
 310   Int_t nOStriplets = arrayOppositeSign->GetEntriesFast();
 311
 312   Int_t nDplusLS=0;
 313   Int_t nLikeSign = arrayLikeSign->GetEntriesFast();
 314   Int_t index;
 315
 316   for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
 317     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
 318     Bool_t unsetvtx=kFALSE;
 319     if(!d->GetOwnPrimaryVtx()) {
 320       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 321       unsetvtx=kTRUE;
 322     }
 323     if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate))nDplusLS++;
 324     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 325   }
 326
 327  Float_t wei2=0;
 328  if(nLikeSign!=0)wei2 = (Float_t)nOStriplets/(Float_t)nLikeSign;
 329  Float_t wei3=0;
 330  if(nDplusLS!=0)wei3 = (Float_t)nDplusOS/(Float_t)nDplusLS;
 331
 332  // loop over like sign candidates
 333   for(Int_t iLikeSign = 0; iLikeSign < nLikeSign; iLikeSign++) {
 334     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)arrayLikeSign->UncheckedAt(iLikeSign);
 335     Bool_t unsetvtx=kFALSE;
 336     if(!d->GetOwnPrimaryVtx()) {
 337       d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 338       unsetvtx=kTRUE;
 339     }
 340
 341     if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate)){
 342
 343       //set tight cuts values
 344       Int_t iPtBin=-1;
 345       Double_t ptCand = d->Pt();
 346       for(Int_t ibin=0;ibin<fNPtBins&&iPtBin<0&&ptCand>fArrayBinLimits[0]&&ptCand<fArrayBinLimits[fNPtBins];ibin++){
 347         if(ptCand<fArrayBinLimits[ibin+1])iPtBin=ibin;
 348       }
 349
 350       if(iPtBin<0){
 351         return;
 352       }
 353
 354       Int_t passTightCuts=fRDCutsAnalysis->IsSelected(d,AliRDHFCuts::kCandidate);
 355
 356       Int_t sign= d->GetCharge();
 357       Float_t wei=1;
 358       Float_t wei4=1;
 359       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
 360
 361         wei=3.*(Float_t)nNegTrks/((Float_t)nPosTrks-2.);
 362         wei4=3.*(Float_t)nspcminus/((Float_t)nspcplus-2.);
 363       }
 364
 365       if(sign<0&&nNegTrks>2&&nspcminus>2){
 366         wei=3.*(Float_t)nPosTrks/((Float_t)nNegTrks-2.);
 367         wei4=3.*(Float_t)nspcplus/((Float_t)nspcminus-2.);
 368
 369       }
 370
 371       Float_t invMass = d->InvMassDplus();
 372       Double_t dlen=d->DecayLength();
 373       Double_t cosp=d->CosPointingAngle();
 374       Double_t sumD02=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 375      Double_t dca=d->GetDCA();
 376      Double_t sigvert=d->GetSigmaVert();
 377      Double_t ptmax=0;
 378      for(Int_t i=0;i<3;i++){
 379        if(d->PtProng(i)>ptmax)ptmax=d->PtProng(i);
 380      }
 381
 382       index=GetLSHistoIndex(iPtBin);
 383       fMassHistLS[index]->Fill(invMass,wei);
 384       fMassHistLS[index+1]->Fill(invMass);
 385       fMassHistLS[index+2]->Fill(invMass,wei2);
 386       fMassHistLS[index+3]->Fill(invMass,wei3);
 387       fMassHistLS[index+4]->Fill(invMass,wei4);
 388
 389       Int_t indexcut=GetHistoIndex(iPtBin);
 390       fCosPHistLS[indexcut]->Fill(cosp);
 391       fDLenHistLS[indexcut]->Fill(dlen);
 392       fSumd02HistLS[indexcut]->Fill(sumD02);
 393       fSigVertHistLS[indexcut]->Fill(sigvert);
 394       fPtMaxHistLS[indexcut]->Fill(ptmax);
 395       fDCAHistLS[indexcut]->Fill(dca);
 396
 397       if(passTightCuts){
 398         fMassHistLSTC[index]->Fill(invMass,wei);
 399         fMassHistLSTC[index+1]->Fill(invMass);
 400         fMassHistLSTC[index+2]->Fill(invMass,wei2);
 401         fMassHistLSTC[index+3]->Fill(invMass,wei3);
 402         fMassHistLSTC[index+4]->Fill(invMass,wei4);
 403       }
 404     }
 405     if(unsetvtx) d->UnsetOwnPrimaryVtx();
 406   }
 407
 408   //printf("------------ N. of positive tracks in Event ----- %d \n", nPosTrks);
 409   //printf("------------ N. of negative tracks in Event ----- %d \n", nNegTrks);
 410
 411   //  printf("LS analysis...done\n");
 412
 413 }
 414
 415
 416 //__________________________________________
 417 void AliAnalysisTaskSEDplus::Init(){
 418   //
 419   // Initialization
 420   //
 421   if(fDebug > 1) printf("AnalysisTaskSEDplus::Init() \n");
 422
 423   //PostData(2,fRDCutsloose);//we should then put those cuts in a tlist if we have more than 1
 424   fListCuts=new TList();
 425   AliRDHFCutsDplustoKpipi *production = new AliRDHFCutsDplustoKpipi();
 426   production=fRDCutsProduction;
 427   AliRDHFCutsDplustoKpipi *analysis = new AliRDHFCutsDplustoKpipi();
 428   analysis=fRDCutsAnalysis;
 429
 430   fListCuts->Add(production);
 431   fListCuts->Add(analysis);
 432   PostData(2,fListCuts);
 433
 434   return;
 435 }
 436
 437 //________________________________________________________________________
 438 void AliAnalysisTaskSEDplus::UserCreateOutputObjects()
 439 {
 440   // Create the output container
 441   //
 442   if(fDebug > 1) printf("AnalysisTaskSEDplus::UserCreateOutputObjects() \n");
 443
 444   // Several histograms are more conveniently managed in a TList
 445   fOutput = new TList();
 446   fOutput->SetOwner();
 447   fOutput->SetName("OutputHistos");
 448
 449   TString hisname;
 450   Int_t index=0;
 451   Int_t indexLS=0;
 452   Int_t nbins=GetNBinsHistos();
 453   for(Int_t i=0;i<fNPtBins;i++){
 454
 455     index=GetHistoIndex(i);
 456     indexLS=GetLSHistoIndex(i);
 457
 458     hisname.Form("hMassPt%d",i);
 459     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 460     fMassHist[index]->Sumw2();
 461     hisname.Form("hCosPAllPt%d",i);
 462     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
 463     fCosPHist[index]->Sumw2();
 464     hisname.Form("hDLenAllPt%d",i);
 465     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
 466     fDLenHist[index]->Sumw2();
 467     hisname.Form("hSumd02AllPt%d",i);
 468     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
 469     fSumd02Hist[index]->Sumw2();
 470     hisname.Form("hSigVertAllPt%d",i);
 471     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 472     fSigVertHist[index]->Sumw2();
 473     hisname.Form("hPtMaxAllPt%d",i);
 474     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 475     fPtMaxHist[index]->Sumw2();
 476
 477     hisname.Form("hDCAAllPt%d",i);
 478     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 479     fDCAHist[index]->Sumw2();
 480
 481
 482
 483     hisname.Form("hMassPt%dTC",i);
 484     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 485     fMassHistTC[index]->Sumw2();
 486     hisname.Form("hMassPt%dTCPlus",i);
 487     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 488     fMassHistTCPlus[index]->Sumw2();
 489     hisname.Form("hMassPt%dTCMinus",i);
 490     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 491     fMassHistTCMinus[index]->Sumw2();
 492
 493
 494
 495
 496     hisname.Form("hCosPAllPt%dLS",i);
 497     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
 498     fCosPHistLS[index]->Sumw2();
 499     hisname.Form("hDLenAllPt%dLS",i);
 500     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
 501     fDLenHistLS[index]->Sumw2();
 502     hisname.Form("hSumd02AllPt%dLS",i);
 503     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
 504     fSumd02HistLS[index]->Sumw2();
 505     hisname.Form("hSigVertAllPt%dLS",i);
 506     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 507     fSigVertHistLS[index]->Sumw2();
 508     hisname.Form("hPtMaxAllPt%dLS",i);
 509     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 510     fPtMaxHistLS[index]->Sumw2();
 511
 512     hisname.Form("hDCAAllPt%dLS",i);
 513     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 514     fDCAHistLS[index]->Sumw2();
 515
 516     hisname.Form("hLSPt%dLC",i);
 517     fMassHistLS[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 518     fMassHistLS[indexLS]->Sumw2();
 519
 520     hisname.Form("hLSPt%dTC",i);
 521     fMassHistLSTC[indexLS] = new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 522     fMassHistLSTC[indexLS]->Sumw2();
 523
 524
 525
 526     index=GetSignalHistoIndex(i);
 527     indexLS++;
 528     hisname.Form("hSigPt%d",i);
 529     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 530     fMassHist[index]->Sumw2();
 531     hisname.Form("hCosPSigPt%d",i);
 532     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
 533     fCosPHist[index]->Sumw2();
 534     hisname.Form("hDLenSigPt%d",i);
 535     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
 536     fDLenHist[index]->Sumw2();
 537     hisname.Form("hSumd02SigPt%d",i);
 538     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
 539     fSumd02Hist[index]->Sumw2();
 540     hisname.Form("hSigVertSigPt%d",i);
 541     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 542     fSigVertHist[index]->Sumw2();
 543     hisname.Form("hPtMaxSigPt%d",i);
 544     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 545     fPtMaxHist[index]->Sumw2();
 546
 547     hisname.Form("hDCASigPt%d",i);
 548     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 549     fDCAHist[index]->Sumw2();
 550
 551
 552     hisname.Form("hSigPt%dTC",i);
 553     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 554     fMassHistTC[index]->Sumw2();
 555     hisname.Form("hSigPt%dTCPlus",i);
 556     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 557     fMassHistTCPlus[index]->Sumw2();
 558     hisname.Form("hSigPt%dTCMinus",i);
 559     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 560     fMassHistTCMinus[index]->Sumw2();
 561
 562     hisname.Form("hLSPt%dLCnw",i);
 563     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 564     fMassHistLS[indexLS]->Sumw2();
 565     hisname.Form("hLSPt%dTCnw",i);
 566     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 567     fMassHistLSTC[indexLS]->Sumw2();
 568
 569
 570
 571     hisname.Form("hCosPSigPt%dLS",i);
 572     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
 573     fCosPHistLS[index]->Sumw2();
 574     hisname.Form("hDLenSigPt%dLS",i);
 575     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
 576     fDLenHistLS[index]->Sumw2();
 577     hisname.Form("hSumd02SigPt%dLS",i);
 578     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
 579     fSumd02HistLS[index]->Sumw2();
 580     hisname.Form("hSigVertSigPt%dLS",i);
 581     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 582     fSigVertHistLS[index]->Sumw2();
 583     hisname.Form("hPtMaxSigPt%dLS",i);
 584     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 585     fPtMaxHistLS[index]->Sumw2();
 586
 587     hisname.Form("hDCASigPt%dLS",i);
 588     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 589     fDCAHistLS[index]->Sumw2();
 590
 591
 592
 593     index=GetBackgroundHistoIndex(i);
 594     indexLS++;
 595     hisname.Form("hBkgPt%d",i);
 596     fMassHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 597     fMassHist[index]->Sumw2();
 598     hisname.Form("hCosPBkgPt%d",i);
 599     fCosPHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
 600     fCosPHist[index]->Sumw2();
 601     hisname.Form("hDLenBkgPt%d",i);
 602     fDLenHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
 603     fDLenHist[index]->Sumw2();
 604     hisname.Form("hSumd02BkgPt%d",i);
 605     fSumd02Hist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
 606     fSumd02Hist[index]->Sumw2();
 607     hisname.Form("hSigVertBkgPt%d",i);
 608     fSigVertHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 609     fSigVertHist[index]->Sumw2();
 610     hisname.Form("hPtMaxBkgPt%d",i);
 611     fPtMaxHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 612     fPtMaxHist[index]->Sumw2();
 613
 614     hisname.Form("hDCABkgPt%d",i);
 615     fDCAHist[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 616     fDCAHist[index]->Sumw2();
 617
 618
 619     hisname.Form("hBkgPt%dTC",i);
 620     fMassHistTC[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 621     fMassHistTC[index]->Sumw2();
 622     hisname.Form("hBkgPt%dTCPlus",i);
 623     fMassHistTCPlus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 624     fMassHistTCPlus[index]->Sumw2();
 625     hisname.Form("hBkgPt%dTCMinus",i);
 626     fMassHistTCMinus[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 627     fMassHistTCMinus[index]->Sumw2();
 628
 629     hisname.Form("hLSPt%dLCntrip",i);
 630     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 631     fMassHistLS[indexLS]->Sumw2();
 632     hisname.Form("hLSPt%dTCntrip",i);
 633     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 634     fMassHistLSTC[indexLS]->Sumw2();
 635
 636
 637     hisname.Form("hCosPBkgPt%dLS",i);
 638     fCosPHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,1.);
 639     fCosPHistLS[index]->Sumw2();
 640     hisname.Form("hDLenBkgPt%dLS",i);
 641     fDLenHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.5);
 642     fDLenHistLS[index]->Sumw2();
 643     hisname.Form("hSumd02BkgPt%dLS",i);
 644     fSumd02HistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,1.);
 645     fSumd02HistLS[index]->Sumw2();
 646     hisname.Form("hSigVertBkgPt%dLS",i);
 647     fSigVertHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 648     fSigVertHistLS[index]->Sumw2();
 649     hisname.Form("hPtMaxBkgPt%dLS",i);
 650     fPtMaxHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.5,5.);
 651     fPtMaxHistLS[index]->Sumw2();
 652     hisname.Form("hDCABkgPt%dLS",i);
 653     fDCAHistLS[index]=new TH1F(hisname.Data(),hisname.Data(),nbins,0.,0.1);
 654     fDCAHistLS[index]->Sumw2();
 655
 656
 657     indexLS++;
 658     hisname.Form("hLSPt%dLCntripsinglecut",i);
 659     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 660     fMassHistLS[indexLS]->Sumw2();
 661     hisname.Form("hLSPt%dTCntripsinglecut",i);
 662     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 663     fMassHistLSTC[indexLS]->Sumw2();
 664
 665     indexLS++;
 666     hisname.Form("hLSPt%dLCspc",i);
 667     fMassHistLS[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 668     fMassHistLS[indexLS]->Sumw2();
 669     hisname.Form("hLSPt%dTCspc",i);
 670     fMassHistLSTC[indexLS]=new TH1F(hisname.Data(),hisname.Data(),nbins,fLowmasslimit,fUpmasslimit);
 671     fMassHistLSTC[indexLS]->Sumw2();
 672   }
 673
 674   for(Int_t i=0; i<3*fNPtBins; i++){
 675     fOutput->Add(fMassHist[i]);
 676      fOutput->Add(fCosPHist[i]);
 677     fOutput->Add(fDLenHist[i]);
 678     fOutput->Add(fSumd02Hist[i]);
 679     fOutput->Add(fSigVertHist[i]);
 680     fOutput->Add(fPtMaxHist[i]);
 681     fOutput->Add(fDCAHist[i]);
 682     fOutput->Add(fMassHistTC[i]);
 683     fOutput->Add(fMassHistTCPlus[i]);
 684     fOutput->Add(fMassHistTCMinus[i]);
 685   }
 686   for(Int_t i=0; i<3*fNPtBins&&fDoLS; i++){
 687     fOutput->Add(fCosPHistLS[i]);
 688     fOutput->Add(fDLenHistLS[i]);
 689     fOutput->Add(fSumd02HistLS[i]);
 690     fOutput->Add(fSigVertHistLS[i]);
 691     fOutput->Add(fPtMaxHistLS[i]);
 692     fOutput->Add(fDCAHistLS[i]);
 693   }
 694   for(Int_t i=0; i<5*fNPtBins&&fDoLS; i++){
 695     fOutput->Add(fMassHistLS[i]);
 696     fOutput->Add(fMassHistLSTC[i]);
 697   }
 698
 699
 700   fHistNEvents = new TH1F("fHistNEvents", "number of events ",6,-0.5,5.5);
 701    fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
 702   fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvents with good vertex");
 703   fHistNEvents->GetXaxis()->SetBinLabel(3,"no. of  pileup events");
 704   fHistNEvents->GetXaxis()->SetBinLabel(4,"no. of candidate");
 705   fHistNEvents->GetXaxis()->SetBinLabel(5,"no. of D+ after loose cuts");
 706   fHistNEvents->GetXaxis()->SetBinLabel(6,"no. of D+ after tight cuts");
 707
 708   fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
 709
 710   fHistNEvents->Sumw2();
 711   fHistNEvents->SetMinimum(0);
 712   fOutput->Add(fHistNEvents);
 713
 714   fPtVsMass=new TH2F("hPtVsMass","PtVsMass (prod. cuts)",nbins,fLowmasslimit,fUpmasslimit,40,0.,20.);
 715   fPtVsMassTC=new TH2F("hPtVsMassTC","PtVsMass (analysis cuts)",nbins,fLowmasslimit,fUpmasslimit,40,0.,20.);
 716   fYVsPt=new TH2F("hYVsPt","YvsPt (prod. cuts)",40,0.,20.,80,-2.,2.);
 717   fYVsPtTC=new TH2F("hYVsPtTC","YvsPt (analysis cuts)",40,0.,20.,80,-2.,2.);
 718   fYVsPtSig=new TH2F("hYVsPtSig","YvsPt (MC, only sig., prod. cuts)",40,0.,20.,80,-2.,2.);
 719   fYVsPtSigTC=new TH2F("hYVsPtSigTC","YvsPt (MC, only Sig, analysis cuts)",40,0.,20.,80,-2.,2.);
 720
 721   fOutput->Add(fPtVsMass);
 722   fOutput->Add(fPtVsMassTC);
 723   fOutput->Add(fYVsPt);
 724   fOutput->Add(fYVsPtTC);
 725   fOutput->Add(fYVsPtSig);
 726   fOutput->Add(fYVsPtSigTC);
 727
 728
 729   //Counter for Normalization
 730   TString normName="NormalizationCounter";
 731   AliAnalysisDataContainer *cont = GetOutputSlot(3)->GetContainer();
 732   if(cont)normName=(TString)cont->GetName();
 733   fCounter = new AliNormalizationCounter(normName.Data());
 734
 735   if(fFillNtuple){
 736     OpenFile(4); // 4 is the slot number of the ntuple
 737
 738     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");
 739
 740   }
 741
 742   return;
 743 }
 744
 745 //________________________________________________________________________
 746 void AliAnalysisTaskSEDplus::UserExec(Option_t */*option*/)
 747 {
 748   // Execute analysis for current event:
 749   // heavy flavor candidates association to MC truth
 750
 751   AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
 752
 753
 754
 755   TClonesArray *array3Prong = 0;
 756   TClonesArray *arrayLikeSign =0;
 757   if(!aod && AODEvent() && IsStandardAOD()) {
 758     // In case there is an AOD handler writing a standard AOD, use the AOD
 759     // event in memory rather than the input (ESD) event.
 760     aod = dynamic_cast<AliAODEvent*> (AODEvent());
 761      // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 762      // have to taken from the AOD event hold by the AliAODExtension
 763     AliAODHandler* aodHandler = (AliAODHandler*)
 764       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 765     if(aodHandler->GetExtensions()) {
 766       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 767       AliAODEvent *aodFromExt = ext->GetAOD();
 768       array3Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 769       arrayLikeSign=(TClonesArray*)aodFromExt->GetList()->FindObject("LikeSign3Prong");
 770     }
 771   } else {
 772     array3Prong=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
 773     arrayLikeSign=(TClonesArray*)aod->GetList()->FindObject("LikeSign3Prong");
 774   }
 775
 776   if(!array3Prong) {
 777     printf("AliAnalysisTaskSEDplus::UserExec: Charm3Prong branch not found!\n");
 778     return;
 779   }
 780   if(!arrayLikeSign) {
 781     printf("AliAnalysisTaskSEDplus::UserExec: LikeSign3Prong branch not found!\n");
 782     return;
 783   }
 784
 785
 786   // fix for temporary bug in ESDfilter
 787   // the AODs with null vertex pointer didn't pass the PhysSel
 788   if(!aod->GetPrimaryVertex()||TMath::Abs(aod->GetMagneticField())<0.001) return;
 789   fCounter->StoreEvent(aod,fReadMC);
 790    fHistNEvents->Fill(0); // count event
 791   // Post the data already here
 792   if(fRDCutsAnalysis->IsEventSelected(aod))fHistNEvents->Fill(1);
 793
 794   if(fRDCutsAnalysis->GetWhyRejection()==1)fHistNEvents->Fill(2);
 795  PostData(1,fOutput);
 796
 797   TClonesArray *arrayMC=0;
 798   AliAODMCHeader *mcHeader=0;
 799
 800   // AOD primary vertex
 801   AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
 802   //    vtx1->Print();
 803    TString primTitle = vtx1->GetTitle();
 804    //if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0)fHistNEvents->Fill(2);
 805
 806   // load MC particles
 807   if(fReadMC){
 808
 809     arrayMC =  (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
 810     if(!arrayMC) {
 811       printf("AliAnalysisTaskSEDplus::UserExec: MC particles branch not found!\n");
 812       //    return;
 813     }
 814
 815   // load MC header
 816     mcHeader =  (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 817     if(!mcHeader) {
 818     printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
 819     return;
 820     }
 821   }
 822
 823   Int_t n3Prong = array3Prong->GetEntriesFast();
 824   if(fDebug>1) printf("Number of D+->Kpipi: %d\n",n3Prong);
 825
 826
 827   Int_t nOS=0;
 828   Int_t index;
 829   Int_t pdgDgDplustoKpipi[3]={321,211,211};
 830   // 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
 831   //Double_t *cutsDplus = new (Double_t*)fRDCuts->GetCuts();
 832   Int_t nSelectedloose=0,nSelectedtight=0;
 833   for (Int_t i3Prong = 0; i3Prong < n3Prong; i3Prong++) {
 834     AliAODRecoDecayHF3Prong *d = (AliAODRecoDecayHF3Prong*)array3Prong->UncheckedAt(i3Prong);
 835     fHistNEvents->Fill(3);
 836
 837     Bool_t unsetvtx=kFALSE;
 838     if(!d->GetOwnPrimaryVtx()){
 839       d->SetOwnPrimaryVtx(vtx1);
 840       unsetvtx=kTRUE;
 841     }
 842
 843     if(fRDCutsProduction->IsSelected(d,AliRDHFCuts::kCandidate)) {
 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);
 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(3,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(4);
 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(5);
 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     if(unsetvtx) d->UnsetOwnPrimaryVtx();
1050   }
1051   fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
1052   fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
1053
1054   //start LS analysis
1055   if(fDoLS && arrayLikeSign) LSAnalysis(array3Prong,arrayLikeSign,aod,vtx1,nOS);
1056
1057   PostData(1,fOutput);
1058   PostData(3,fCounter);
1059   return;
1060 }
1061
1062
1063
1064 //________________________________________________________________________
1065 void AliAnalysisTaskSEDplus::Terminate(Option_t */*option*/)
1066 {
1067   // Terminate analysis
1068   //
1069   if(fDebug > 1) printf("AnalysisTaskSEDplus: Terminate() \n");
1070
1071   fOutput = dynamic_cast<TList*> (GetOutputData(1));
1072   if (!fOutput) {
1073     printf("ERROR: fOutput not available\n");
1074     return;
1075   }
1076   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
1077   fYVsPt = dynamic_cast<TH2F*>(fOutput->FindObject("hYVsPt"));
1078   fYVsPtTC = dynamic_cast<TH2F*>(fOutput->FindObject("hYVsPtTC"));
1079   fYVsPtSig = dynamic_cast<TH2F*>(fOutput->FindObject("hYVsPtSig"));
1080   fYVsPtSigTC = dynamic_cast<TH2F*>(fOutput->FindObject("hYVsPtSigTC"));
1081   fPtVsMass = dynamic_cast<TH2F*>(fOutput->FindObject("hPtVsMass"));
1082   fPtVsMassTC = dynamic_cast<TH2F*>(fOutput->FindObject("hPtVsMassTC"));
1083
1084   TString hisname;
1085   Int_t index=0;
1086
1087
1088   Int_t indexLS=0;
1089   for(Int_t i=0;i<fNPtBins;i++){
1090     index=GetHistoIndex(i);
1091     if(fDoLS)indexLS=GetLSHistoIndex(i);
1092     hisname.Form("hMassPt%d",i);
1093     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1094     hisname.Form("hCosPAllPt%d",i);
1095     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1096     hisname.Form("hDLenAllPt%d",i);
1097     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1098     hisname.Form("hSumd02AllPt%d",i);
1099     fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1100     hisname.Form("hSigVertAllPt%d",i);
1101     fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1102     hisname.Form("hPtMaxAllPt%d",i);
1103     fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1104     hisname.Form("hDCAAllPt%d",i);
1105     fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1106     hisname.Form("hMassPt%dTC",i);
1107     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1108     hisname.Form("hMassPt%dTCPlus",i);
1109     fMassHistTCPlus[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1110     hisname.Form("hMassPt%dTCMinus",i);
1111     fMassHistTCMinus[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1112     if(fDoLS){
1113       hisname.Form("hLSPt%dLC",i);
1114       fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1115       hisname.Form("hCosPAllPt%dLS",i);
1116       fCosPHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1117       hisname.Form("hDLenAllPt%dLS",i);
1118       fDLenHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1119       hisname.Form("hSumd02AllPt%dLS",i);
1120       fSumd02HistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1121       hisname.Form("hSigVertAllPt%dLS",i);
1122       fSigVertHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1123       hisname.Form("hPtMaxAllPt%dLS",i);
1124       fPtMaxHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1125       hisname.Form("hDCAAllPt%dLS",i);
1126       fDCAHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1127
1128       hisname.Form("hLSPt%dTC",i);
1129       fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1130
1131     }
1132
1133     index=GetSignalHistoIndex(i);
1134     if(fDoLS)indexLS++;
1135     hisname.Form("hSigPt%d",i);
1136     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1137     hisname.Form("hCosPSigPt%d",i);
1138     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1139     hisname.Form("hDLenSigPt%d",i);
1140     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1141     hisname.Form("hSumd02SigPt%d",i);
1142     fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1143     hisname.Form("hSigVertSigPt%d",i);
1144     fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1145     hisname.Form("hPtMaxSigPt%d",i);
1146     fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1147     hisname.Form("hDCASigPt%d",i);
1148     fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1149
1150     hisname.Form("hSigPt%dTC",i);
1151     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1152     hisname.Form("hSigPt%dTCPlus",i);
1153     fMassHistTCPlus[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1154     hisname.Form("hSigPt%dTCMinus",i);
1155     fMassHistTCMinus[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1156     if(fDoLS){
1157       hisname.Form("hLSPt%dLCnw",i);
1158       fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1159       hisname.Form("hCosPSigPt%dLS",i);
1160       fCosPHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1161       hisname.Form("hDLenSigPt%dLS",i);
1162       fDLenHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1163       hisname.Form("hSumd02SigPt%dLS",i);
1164       fSumd02HistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1165       hisname.Form("hSigVertSigPt%dLS",i);
1166       fSigVertHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1167       hisname.Form("hPtMaxSigPt%dLS",i);
1168       fPtMaxHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1169       hisname.Form("hDCASigPt%dLS",i);
1170       fDCAHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1171
1172       hisname.Form("hLSPt%dTCnw",i);
1173       fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1174     }
1175
1176     index=GetBackgroundHistoIndex(i);
1177     if(fDoLS)indexLS++;
1178     hisname.Form("hBkgPt%d",i);
1179     fMassHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1180     hisname.Form("hCosPBkgPt%d",i);
1181     fCosPHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1182     hisname.Form("hDLenBkgPt%d",i);
1183     fDLenHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1184     hisname.Form("hSumd02BkgPt%d",i);
1185     fSumd02Hist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1186     hisname.Form("hSigVertBkgPt%d",i);
1187     fSigVertHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1188     hisname.Form("hPtMaxBkgPt%d",i);
1189     fPtMaxHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1190     hisname.Form("hDCABkgPt%d",i);
1191     fDCAHist[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1192     hisname.Form("hBkgPt%dTC",i);
1193     fMassHistTC[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1194     hisname.Form("hBkgPt%dTCPlus",i);
1195     fMassHistTCPlus[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1196     hisname.Form("hBkgPt%dTCMinus",i);
1197     fMassHistTCMinus[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1198     if(fDoLS){
1199       hisname.Form("hLSPt%dLCntrip",i);
1200       fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1201
1202
1203       hisname.Form("hCosPBkgPt%dLS",i);
1204       fCosPHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1205       hisname.Form("hDLenBkgPt%dLS",i);
1206       fDLenHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1207       hisname.Form("hSumd02BkgPt%dLS",i);
1208       fSumd02HistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1209       hisname.Form("hSigVertBkgPt%dLS",i);
1210       fSigVertHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1211       hisname.Form("hPtMaxBkgPt%dLS",i);
1212       fPtMaxHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1213       hisname.Form("hDCABkgPt%dLS",i);
1214       fDCAHistLS[index]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1215
1216       hisname.Form("hLSPt%dTCntrip",i);
1217       fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1218
1219
1220       indexLS++;
1221       hisname.Form("hLSPt%dLCntripsinglecut",i);
1222       fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1223
1224       hisname.Form("hLSPt%dTCntripsinglecut",i);
1225       fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1226
1227
1228       indexLS++;
1229       hisname.Form("hLSPt%dLCspc",i);
1230       fMassHistLS[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1231
1232       hisname.Form("hLSPt%dTCspc",i);
1233       fMassHistLSTC[indexLS]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1234     }
1235
1236   }
1237   fCounter = dynamic_cast<AliNormalizationCounter*>(GetOutputData(3));
1238
1239   if(fFillNtuple){
1240     fNtupleDplus = dynamic_cast<TNtuple*>(GetOutputData(4));
1241   }
1242
1243   TCanvas *c1=new TCanvas("c1","D+ invariant mass distribution",500,500);
1244   c1->cd();
1245   TH1F *hMassPt=(TH1F*)fOutput->FindObject("hMassPt3TC");
1246   hMassPt->SetLineColor(kBlue);
1247   hMassPt->SetXTitle("M[GeV/c^{2}]");
1248   hMassPt->Draw();
1249
1250   return;
1251 }