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