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