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