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