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