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