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