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