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