]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG3/vertexingHF/AliAnalysisTaskSESelectHF4Prong.cxx
Update (Fabio)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisTaskSESelectHF4Prong.cxx
index 6e8d7d9530fb7452bfedcb311c3eea1b64967c9a..c6aa305c520422494740cc419d141bdf71c64f48 100644 (file)
@@ -59,6 +59,7 @@ AliAnalysisTaskSE(),
 fVerticesHFTClArr(0),
 fCharm4ProngTClArr(0),
 fSelected(0),
+fMCTruth(0),
 fOutput(0),
 fOutput2(0),
 fOutput3(0),
@@ -95,11 +96,21 @@ fhInvMassSumAll10MevBin5(0),
 fhInvMassD0Sum5MevBin5(0),
 fhInvMassD0barSum5MevBin5(0),
 fhInvMassSumAll5MevBin5(0),
-fhInvMassMultipleOnlyBin1(0),
-fhInvMassMultipleOnlyBin2(0),
-fhInvMassMultipleOnlyBin3(0),
-fhInvMassMultipleOnlyBin4(0),
-fhInvMassMultipleOnlyBin5(0),
+fhReflBin1(0),
+fhReflBin2(0),
+fhReflBin3(0),
+fhReflBin4(0),
+fhReflBin5(0),
+fhReflD0Bin1(0),
+fhReflD0Bin2(0),
+fhReflD0Bin3(0),
+fhReflD0Bin4(0),
+fhReflD0Bin5(0),
+fhReflD0barBin1(0),
+fhReflD0barBin2(0),
+fhReflD0barBin3(0),
+fhReflD0barBin4(0),
+fhReflD0barBin5(0),
 fScatterP4PID(0),
 fPtVsY(0),
 fPtVsYAll(0),
@@ -128,12 +139,13 @@ fCuts(0)
   // Default constructor
 }
 
-//________________________________________________________________________
+//_______________________________________________________________________
 AliAnalysisTaskSESelectHF4Prong::AliAnalysisTaskSESelectHF4Prong(const char *name,AliRDHFCutsD0toKpipipi* cuts):
 AliAnalysisTaskSE(name),
 fVerticesHFTClArr(0),
 fCharm4ProngTClArr(0),
 fSelected(0),
+fMCTruth(0),
 fOutput(0),
 fOutput2(0),
 fOutput3(0),
@@ -170,11 +182,21 @@ fhInvMassSumAll10MevBin5(0),
 fhInvMassD0Sum5MevBin5(0),
 fhInvMassD0barSum5MevBin5(0),
 fhInvMassSumAll5MevBin5(0),
-fhInvMassMultipleOnlyBin1(0),
-fhInvMassMultipleOnlyBin2(0),
-fhInvMassMultipleOnlyBin3(0),
-fhInvMassMultipleOnlyBin4(0),
-fhInvMassMultipleOnlyBin5(0),
+fhReflBin1(0),
+fhReflBin2(0),
+fhReflBin3(0),
+fhReflBin4(0),
+fhReflBin5(0),
+fhReflD0Bin1(0),
+fhReflD0Bin2(0),
+fhReflD0Bin3(0),
+fhReflD0Bin4(0),
+fhReflD0Bin5(0),
+fhReflD0barBin1(0),
+fhReflD0barBin2(0),
+fhReflD0barBin3(0),
+fhReflD0barBin4(0),
+fhReflD0barBin5(0),
 fScatterP4PID(0),
 fPtVsY(0),
 fPtVsYAll(0),
@@ -261,6 +283,7 @@ void AliAnalysisTaskSESelectHF4Prong::Init()
   // Initialization
 
   if(fDebug > 1) printf("AnalysisTaskSESelectHF4Prong::Init() \n");
+  PrintPtBinHandMCFlag();
 
   return;
 }
@@ -272,6 +295,8 @@ void AliAnalysisTaskSESelectHF4Prong::UserCreateOutputObjects()
   //
   if(fDebug > 1) printf("AnalysisTaskSESelectHF4Prong::UserCreateOutputObjects() \n");
 
+  if(fDebug > 1) PrintPtBinHandMCFlag();
+
   fVerticesHFTClArr = new TClonesArray("AliAODVertex", 0);
   fVerticesHFTClArr->SetName("VerticesHF");
   AddAODBranch("TClonesArray", &fVerticesHFTClArr);
@@ -448,30 +473,80 @@ void AliAnalysisTaskSESelectHF4Prong::UserCreateOutputObjects()
   fhInvMassSumAll5MevBin5->SetMinimum(0);
   fOutput5->Add(fhInvMassSumAll5MevBin5);
 
-  fhInvMassMultipleOnlyBin1 = new TH1F("fhInvMassMultipleOnlyBin1", "D0/D0bar invariant mass Bin1 (good hyps) - Multple hyps accepted only; Inv. mass [GeV]; Entries/10 MeV",120,1.6,2.2);
-  fhInvMassMultipleOnlyBin1->Sumw2(); //Create structure to store sum of squares of weights
-  fhInvMassMultipleOnlyBin1->SetMinimum(0);
-  fOutput->Add(fhInvMassMultipleOnlyBin1);
-
-  fhInvMassMultipleOnlyBin2 = new TH1F("fhInvMassMultipleOnlyBin2", "D0/D0bar invariant mass Bin2 (good hyps) - Multple hyps accepted only; Inv. mass [GeV]; Entries/10 MeV",120,1.6,2.2);
-  fhInvMassMultipleOnlyBin2->Sumw2(); //Create structure to store sum of squares of weights
-  fhInvMassMultipleOnlyBin2->SetMinimum(0);
-  fOutput2->Add(fhInvMassMultipleOnlyBin2);
-
-  fhInvMassMultipleOnlyBin3 = new TH1F("fhInvMassMultipleOnlyBin3", "D0/D0bar invariant mass Bin3 (good hyps) - Multple hyps accepted only; Inv. mass [GeV]; Entries/10 MeV",120,1.6,2.2);
-  fhInvMassMultipleOnlyBin3->Sumw2(); //Create structure to store sum of squares of weights
-  fhInvMassMultipleOnlyBin3->SetMinimum(0);
-  fOutput3->Add(fhInvMassMultipleOnlyBin3);
-
-  fhInvMassMultipleOnlyBin4 = new TH1F("fhInvMassMultipleOnlyBin4", "D0/D0bar invariant mass Bin4 (good hyps) - Multple hyps accepted only; Inv. mass [GeV]; Entries/10 MeV",120,1.6,2.2);
-  fhInvMassMultipleOnlyBin4->Sumw2(); //Create structure to store sum of squares of weights
-  fhInvMassMultipleOnlyBin4->SetMinimum(0);
-  fOutput4->Add(fhInvMassMultipleOnlyBin4);
-
-  fhInvMassMultipleOnlyBin5 = new TH1F("fhInvMassMultipleOnlyBin5", "D0/D0bar invariant mass Bin5 (good hyps) - Multple hyps accepted only; Inv. mass [GeV]; Entries/10 MeV",120,1.6,2.2);
-  fhInvMassMultipleOnlyBin5->Sumw2(); //Create structure to store sum of squares of weights
-  fhInvMassMultipleOnlyBin5->SetMinimum(0);
-  fOutput5->Add(fhInvMassMultipleOnlyBin5);
+  fhReflBin1 = new TH2F("fhReflBin1", "Invariant Mass Histogram for reflections; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,8,0.,8.);
+  fhReflBin1->Sumw2(); //Create structure to store sum of squares of weights
+  fhReflBin1->SetMinimum(0);
+  fOutput->Add(fhReflBin1);
+
+  fhReflBin2 = new TH2F("fhReflBin2", "Invariant Mass Histogram for reflections; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,8,0.,8.);
+  fhReflBin2->Sumw2(); //Create structure to store sum of squares of weights
+  fhReflBin2->SetMinimum(0);
+  fOutput2->Add(fhReflBin2);
+
+  fhReflBin3 = new TH2F("fhReflBin3", "Invariant Mass Histogram for reflections; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,8,0.,8.);
+  fhReflBin3->Sumw2(); //Create structure to store sum of squares of weights
+  fhReflBin3->SetMinimum(0);
+  fOutput3->Add(fhReflBin3);
+
+  fhReflBin4 = new TH2F("fhReflBin4", "Invariant Mass Histogram for reflections; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,8,0.,8.);
+  fhReflBin4->Sumw2(); //Create structure to store sum of squares of weights
+  fhReflBin4->SetMinimum(0);
+  fOutput4->Add(fhReflBin4);
+
+  fhReflBin5 = new TH2F("fhReflBin5", "Invariant Mass Histogram for reflections; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,8,0.,8.);
+  fhReflBin5->Sumw2(); //Create structure to store sum of squares of weights
+  fhReflBin5->SetMinimum(0);
+  fOutput5->Add(fhReflBin5);
+
+  fhReflD0Bin1 = new TH2F("fhReflD0Bin1", "Invariant Mass Histogram for reflections - D0 hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
+  fhReflD0Bin1->Sumw2(); //Create structure to store sum of squares of weights
+  fhReflD0Bin1->SetMinimum(0);
+  fOutput->Add(fhReflD0Bin1);
+
+  fhReflD0Bin2 = new TH2F("fhReflD0Bin2", "Invariant Mass Histogram for reflections - D0 hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
+  fhReflD0Bin2->Sumw2(); //Create structure to store sum of squares of weights
+  fhReflD0Bin2->SetMinimum(0);
+  fOutput2->Add(fhReflD0Bin2);
+
+  fhReflD0Bin3 = new TH2F("fhReflD0Bin3", "Invariant Mass Histogram for reflections - D0 hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
+  fhReflD0Bin3->Sumw2(); //Create structure to store sum of squares of weights
+  fhReflD0Bin3->SetMinimum(0);
+  fOutput3->Add(fhReflD0Bin3);
+
+  fhReflD0Bin4 = new TH2F("fhReflD0Bin4", "Invariant Mass Histogram for reflections - D0 hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
+  fhReflD0Bin4->Sumw2(); //Create structure to store sum of squares of weights
+  fhReflD0Bin4->SetMinimum(0);
+  fOutput4->Add(fhReflD0Bin4);
+
+  fhReflD0Bin5 = new TH2F("fhReflD0Bin5", "Invariant Mass Histogram for reflections - D0 hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
+  fhReflD0Bin5->Sumw2(); //Create structure to store sum of squares of weights
+  fhReflD0Bin5->SetMinimum(0);
+  fOutput5->Add(fhReflD0Bin5);
+
+  fhReflD0barBin1 = new TH2F("fhReflD0barBin1", "Invariant Mass Histogram for reflections - D0bar hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
+  fhReflD0barBin1->Sumw2(); //Create structure to store sum of squares of weights
+  fhReflD0barBin1->SetMinimum(0);
+  fOutput->Add(fhReflD0barBin1);
+
+  fhReflD0barBin2 = new TH2F("fhReflD0barBin2", "Invariant Mass Histogram for reflections - D0bar hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
+  fhReflD0barBin2->Sumw2(); //Create structure to store sum of squares of weights
+  fhReflD0barBin2->SetMinimum(0);
+  fOutput2->Add(fhReflD0barBin2);
+
+  fhReflD0barBin3 = new TH2F("fhReflD0barBin3", "Invariant Mass Histogram for reflections - D0bar hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
+  fhReflD0barBin3->Sumw2(); //Create structure to store sum of squares of weights
+  fhReflD0barBin3->SetMinimum(0);
+  fOutput3->Add(fhReflD0barBin3);
+
+  fhReflD0barBin4 = new TH2F("fhReflD0barBin4", "Invariant Mass Histogram for reflections - D0bar hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
+  fhReflD0barBin4->Sumw2(); //Create structure to store sum of squares of weights
+  fhReflD0barBin4->SetMinimum(0);
+  fOutput4->Add(fhReflD0barBin4);
+
+  fhReflD0barBin5 = new TH2F("fhReflD0barBin5", "Invariant Mass Histogram for reflections - D0bar hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
+  fhReflD0barBin5->Sumw2(); //Create structure to store sum of squares of weights
+  fhReflD0barBin5->SetMinimum(0);
+  fOutput5->Add(fhReflD0barBin5);
 
   fScatterP4PID = new TH2F("fScatterP4PID", "Transverse momentum of K vs l-s Pi (D0 + D0bar); Pt of K [GeV/c]; Pt of Pi [GeV/c]",500,0.,5.,500,0.,5.);
   fScatterP4PID->SetMinimum(0);
@@ -565,7 +640,7 @@ void AliAnalysisTaskSESelectHF4Prong::UserCreateOutputObjects()
   fPIDSelBin5->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
   fPIDSelBin5->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
 
-  fMultipleHyps = new TH1F("fMultipleHyps", "N. of hyp. accepted for each candidate (accounted N. times)",8,0.,8.);
+  fMultipleHyps = new TH2F("fMultipleHyps", "N. of hyp. accepted for each candidate (accounted N. times)",8,0.,8.,5,0.,5.);
   fMultipleHyps->SetMinimum(0);
   fMultipleHyps->GetXaxis()->SetBinLabel(1,"1 (PPR)");
   fMultipleHyps->GetXaxis()->SetBinLabel(2,"2 (PPR)");
@@ -575,8 +650,13 @@ void AliAnalysisTaskSESelectHF4Prong::UserCreateOutputObjects()
   fMultipleHyps->GetXaxis()->SetBinLabel(6,"2 (PPR+PID)");
   fMultipleHyps->GetXaxis()->SetBinLabel(7,"3 (PPR+PID)");
   fMultipleHyps->GetXaxis()->SetBinLabel(8,"4 (PPR+PID)");
+  fMultipleHyps->GetYaxis()->SetBinLabel(1,"PtBin 1");
+  fMultipleHyps->GetYaxis()->SetBinLabel(2,"PtBin 2");
+  fMultipleHyps->GetYaxis()->SetBinLabel(3,"PtBin 3");
+  fMultipleHyps->GetYaxis()->SetBinLabel(4,"PtBin 4");
+  fMultipleHyps->GetYaxis()->SetBinLabel(5,"PtBin 5");
 
-  fMultipleHypsType = new TH1F("fMultipleHypsType", "Type of hyp. accepted for each candidate",8,0.,8.);
+  fMultipleHypsType = new TH2F("fMultipleHypsType", "Type of hyp. accepted for each candidate",8,0.,8.,5,0.,5.);
   fMultipleHypsType->SetMinimum(0);
   fMultipleHypsType->GetXaxis()->SetBinLabel(1,"D0");
   fMultipleHypsType->GetXaxis()->SetBinLabel(2,"D0bar");
@@ -586,6 +666,11 @@ void AliAnalysisTaskSESelectHF4Prong::UserCreateOutputObjects()
   fMultipleHypsType->GetXaxis()->SetBinLabel(6,"2D0+D0bar");
   fMultipleHypsType->GetXaxis()->SetBinLabel(7,"D0+2D0bar");
   fMultipleHypsType->GetXaxis()->SetBinLabel(8,"2D0+2D0bar");
+  fMultipleHypsType->GetYaxis()->SetBinLabel(1,"PtBin 1");
+  fMultipleHypsType->GetYaxis()->SetBinLabel(2,"PtBin 2");
+  fMultipleHypsType->GetYaxis()->SetBinLabel(3,"PtBin 3");
+  fMultipleHypsType->GetYaxis()->SetBinLabel(4,"PtBin 4");
+  fMultipleHypsType->GetYaxis()->SetBinLabel(5,"PtBin 5");
 
   fOutputC->Add(fMultipleHyps);
   fOutputC->Add(fMultipleHypsType);
@@ -692,8 +777,6 @@ void AliAnalysisTaskSESelectHF4Prong::UserExec(Option_t */*option*/)
    Double_t dcaquad1 = dIn->GetDCA(2);
    Double_t dcaquad2 = dIn->GetDCA(5);
 
-   Double_t ptBinH[6] = {2.5,4.5,6.0,8.0,12.0,25.0};  //bin i has pt between values i and i+1
-
    fCutDCA->Fill(dca);
    fCutDCA3->Fill(dcatrip);
    fCutDCA2->Fill(dcaquad1);
@@ -728,14 +811,13 @@ void AliAnalysisTaskSESelectHF4Prong::UserExec(Option_t */*option*/)
    if(fSelected==1||fSelected==2||fSelected==3) {
       fPtVsY->Fill(ptPart,yPart);
       fPIDSel->Fill(0);
-        if (ptPart >= ptBinH[0] && ptPart < ptBinH[1])             fPIDSelBin1->Fill(0);
-        else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fPIDSelBin2->Fill(0);
-        else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fPIDSelBin3->Fill(0);
-        else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fPIDSelBin4->Fill(0);
-        else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fPIDSelBin5->Fill(0);
+        if (ptPart >= fPtBinH[0] && ptPart < fPtBinH[1])      fPIDSelBin1->Fill(0);
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fPIDSelBin2->Fill(0);
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fPIDSelBin3->Fill(0);
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fPIDSelBin4->Fill(0);
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fPIDSelBin5->Fill(0);
         }
       
-
    PostData(6,fOutputC);
    
    //selection
@@ -760,6 +842,7 @@ void AliAnalysisTaskSESelectHF4Prong::UserExec(Option_t */*option*/)
       dIn->InvMassD0bar(fmassD0bar);
 
       //fill histos (combining selection form cuts & PID (with rho information))
+      Double_t binPt = -1;
       Int_t hypD01 = 0, hypD02 = 0, hypD0bar1 = 0, hypD0bar2 = 0;
       Int_t pid1 = 0, pid2 = 0, pidbar1 = 0, pidbar2 = 0;
       Int_t pidSelection = fCuts->IsSelectedFromPID(dIn, &pid1, &pid2, &pidbar1, &pidbar2);
@@ -767,18 +850,18 @@ void AliAnalysisTaskSESelectHF4Prong::UserExec(Option_t */*option*/)
    //number of CANDIDATES (regardless of hypotheses) passing PPR + PID - PAY ATTENTION: hypoth. for PID and PPR may not be the same!
    if (pidSelection > 0) {
         fPIDSel->Fill(1);
-        if (ptPart >= ptBinH[0] && ptPart < ptBinH[1])             fPIDSelBin1->Fill(1);
-        else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fPIDSelBin2->Fill(1);
-        else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fPIDSelBin3->Fill(1);
-        else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fPIDSelBin4->Fill(1);
-        else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fPIDSelBin5->Fill(1);
+        if (ptPart >= fPtBinH[0] && ptPart < fPtBinH[1])      {fPIDSelBin1->Fill(1); binPt = 0;}
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) {fPIDSelBin2->Fill(1); binPt = 1;}
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) {fPIDSelBin3->Fill(1); binPt = 2;}
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) {fPIDSelBin4->Fill(1); binPt = 3;}
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) {fPIDSelBin5->Fill(1); binPt = 4;}
       }
    
           //number of hypoteses accepted per candidate after PPR
-          if(selD01+selD02+selD0bar1+selD0bar2 == 1) fMultipleHyps->Fill(0);
-          if(selD01+selD02+selD0bar1+selD0bar2 == 2) fMultipleHyps->Fill(1);
-          if(selD01+selD02+selD0bar1+selD0bar2 == 3) fMultipleHyps->Fill(2);
-          if(selD01+selD02+selD0bar1+selD0bar2 == 4) fMultipleHyps->Fill(3);
+          if(selD01+selD02+selD0bar1+selD0bar2 == 1) fMultipleHyps->Fill((Double_t)0,binPt);
+          if(selD01+selD02+selD0bar1+selD0bar2 == 2) fMultipleHyps->Fill((Double_t)1,binPt);
+          if(selD01+selD02+selD0bar1+selD0bar2 == 3) fMultipleHyps->Fill((Double_t)2,binPt);
+          if(selD01+selD02+selD0bar1+selD0bar2 == 4) fMultipleHyps->Fill((Double_t)3,binPt);
 
    //combine PPD + PID cuts
    if (selD01 == 1 && pid1==1) hypD01 = 1;
@@ -789,101 +872,56 @@ void AliAnalysisTaskSESelectHF4Prong::UserExec(Option_t */*option*/)
    //number of CANDIDATES (regardless of hypotheses) passing PPR + PID - PAY ATTENTION: hypoth. for PID and PPR must match (at least one)!
    if (hypD01 == 1 || hypD02 == 1 || hypD0bar1 == 1 || hypD0bar2 == 1) {
         fPIDSel->Fill(2);
-        if (ptPart >= ptBinH[0] && ptPart < ptBinH[1])             fPIDSelBin1->Fill(2);
-        else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fPIDSelBin2->Fill(2);
-        else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fPIDSelBin3->Fill(2);
-        else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fPIDSelBin4->Fill(2);
-        else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fPIDSelBin5->Fill(2);
+        if (ptPart >= fPtBinH[0] && ptPart < fPtBinH[1])      {fPIDSelBin1->Fill(2); binPt = 0;}
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) {fPIDSelBin2->Fill(2); binPt = 1;}
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) {fPIDSelBin3->Fill(2); binPt = 2;}
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) {fPIDSelBin4->Fill(2); binPt = 3;}
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) {fPIDSelBin5->Fill(2); binPt = 4;}
       }
 
           //number of hypoteses accepted per candidate after PPR and PID
-          if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 1) fMultipleHyps->Fill(4);
-          if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 2) fMultipleHyps->Fill(5);
-          if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 3) fMultipleHyps->Fill(6);
-          if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 4) fMultipleHyps->Fill(7);
+          if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 1) fMultipleHyps->Fill((Double_t)4,binPt);
+          if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 2) fMultipleHyps->Fill((Double_t)5,binPt);
+          if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 3) fMultipleHyps->Fill((Double_t)6,binPt);
+          if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 4) fMultipleHyps->Fill((Double_t)7,binPt);
 
           //type of hypoteses accepted per candidate after PPR and PID
-          if(hypD01+hypD02 == 1 && hypD0bar1+hypD0bar2 == 0) fMultipleHypsType->Fill(0);
-          if(hypD01+hypD02 == 0 && hypD0bar1+hypD0bar2 == 1) fMultipleHypsType->Fill(1);
-          if(hypD01+hypD02 == 2 && hypD0bar1+hypD0bar2 == 0) fMultipleHypsType->Fill(2);
-          if(hypD01+hypD02 == 0 && hypD0bar1+hypD0bar2 == 2) fMultipleHypsType->Fill(3);
-          if(hypD01+hypD02 == 1 && hypD0bar1+hypD0bar2 == 1) fMultipleHypsType->Fill(4);
-          if(hypD01+hypD02 == 2 && hypD0bar1+hypD0bar2 == 1) fMultipleHypsType->Fill(5);
-          if(hypD01+hypD02 == 1 && hypD0bar1+hypD0bar2 == 2) fMultipleHypsType->Fill(6);
-          if(hypD01+hypD02 == 2 && hypD0bar1+hypD0bar2 == 2) fMultipleHypsType->Fill(7);
-
-
-//
-// This section is auxiliary: (multiple hypotheses histos)
-//
- if (hypD01+hypD02+hypD0bar1+hypD0bar2 == 2 || hypD01+hypD02+hypD0bar1+hypD0bar2 == 3 || hypD01+hypD02+hypD0bar1+hypD0bar2 == 4) 
- {
-    
-   if (ptPart > ptBinH[0]) {
-      // D01 hyp.
-      if(hypD01==1) {
-        if (ptPart < ptBinH[1])                            fhInvMassMultipleOnlyBin1->Fill(fmassD0[0]);
-        else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fhInvMassMultipleOnlyBin2->Fill(fmassD0[0]);
-        else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fhInvMassMultipleOnlyBin3->Fill(fmassD0[0]);
-        else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fhInvMassMultipleOnlyBin4->Fill(fmassD0[0]);
-        else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fhInvMassMultipleOnlyBin5->Fill(fmassD0[0]);
-       }
-      // D02 hyp.
-      if(hypD02==1) {
-        if (ptPart < ptBinH[1])                            fhInvMassMultipleOnlyBin1->Fill(fmassD0[1]);
-        else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fhInvMassMultipleOnlyBin2->Fill(fmassD0[1]);
-        else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fhInvMassMultipleOnlyBin3->Fill(fmassD0[1]);
-        else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fhInvMassMultipleOnlyBin4->Fill(fmassD0[1]);
-        else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fhInvMassMultipleOnlyBin5->Fill(fmassD0[1]);
-       }
-      // D0bar1 hyp.
-      if(hypD0bar1==1) {
-       if (ptPart < ptBinH[1])                             fhInvMassMultipleOnlyBin1->Fill(fmassD0bar[0]);
-        else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fhInvMassMultipleOnlyBin2->Fill(fmassD0bar[0]);
-        else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fhInvMassMultipleOnlyBin3->Fill(fmassD0bar[0]);
-        else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fhInvMassMultipleOnlyBin4->Fill(fmassD0bar[0]);
-        else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fhInvMassMultipleOnlyBin5->Fill(fmassD0bar[0]); 
-       } 
-      // D0bar2 hyp.
-      if(hypD0bar2==1) {
-        if (ptPart < ptBinH[1])                            fhInvMassMultipleOnlyBin1->Fill(fmassD0bar[1]);
-        else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fhInvMassMultipleOnlyBin2->Fill(fmassD0bar[1]);
-        else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fhInvMassMultipleOnlyBin3->Fill(fmassD0bar[1]);
-        else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fhInvMassMultipleOnlyBin4->Fill(fmassD0bar[1]);
-        else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fhInvMassMultipleOnlyBin5->Fill(fmassD0bar[1]);
-       }
-     }
-
-   }
-// 
-//end of auxiliary section
-//
-
+          if(hypD01+hypD02 == 1 && hypD0bar1+hypD0bar2 == 0) fMultipleHypsType->Fill((Double_t)0,binPt);
+          if(hypD01+hypD02 == 0 && hypD0bar1+hypD0bar2 == 1) fMultipleHypsType->Fill((Double_t)1,binPt);
+          if(hypD01+hypD02 == 2 && hypD0bar1+hypD0bar2 == 0) fMultipleHypsType->Fill((Double_t)2,binPt);
+          if(hypD01+hypD02 == 0 && hypD0bar1+hypD0bar2 == 2) fMultipleHypsType->Fill((Double_t)3,binPt);
+          if(hypD01+hypD02 == 1 && hypD0bar1+hypD0bar2 == 1) fMultipleHypsType->Fill((Double_t)4,binPt);
+          if(hypD01+hypD02 == 2 && hypD0bar1+hypD0bar2 == 1) fMultipleHypsType->Fill((Double_t)5,binPt);
+          if(hypD01+hypD02 == 1 && hypD0bar1+hypD0bar2 == 2) fMultipleHypsType->Fill((Double_t)6,binPt);
+          if(hypD01+hypD02 == 2 && hypD0bar1+hypD0bar2 == 2) fMultipleHypsType->Fill((Double_t)7,binPt);
+
+   //Call function for reflection analysis
+   if (hypD01+hypD02+hypD0bar1+hypD0bar2 > 0) AnalysisReflection(aodIn, dIn, hypD01, hypD02, hypD0bar1, hypD0bar2);
 
    //All histos are filled if Pt of candidate is greater than minimum of first bin (in this way: bin1+bin2+...binN = whole)
-   if (ptPart > ptBinH[0]) {
+   if (ptPart > fPtBinH[0]) {
     
       // D01 hyp.
       if(hypD01==1) {
        fPtSel->Fill(ptPart);
         fScatterP4PID->Fill(trk1->Pt(),trk3->Pt());
-        if (ptPart < ptBinH[1])                         {fhInvMassD0Sum10MevBin1->Fill(fmassD0[0]);
+        if (ptPart < fPtBinH[1])                         {fhInvMassD0Sum10MevBin1->Fill(fmassD0[0]);
                                                          fhInvMassD0Sum5MevBin1->Fill(fmassD0[0]);
                                                          fhInvMassSumAll10MevBin1->Fill(fmassD0[0]);
                                                          fhInvMassSumAll5MevBin1->Fill(fmassD0[0]);}
-        else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) {fhInvMassD0Sum10MevBin2->Fill(fmassD0[0]);
-                                                            fhInvMassD0Sum5MevBin2->Fill(fmassD0[0]);
-                                                            fhInvMassSumAll10MevBin2->Fill(fmassD0[0]);
-                                                            fhInvMassSumAll5MevBin2->Fill(fmassD0[0]);}
-        else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) {fhInvMassD0Sum10MevBin3->Fill(fmassD0[0]);
-                                                            fhInvMassD0Sum5MevBin3->Fill(fmassD0[0]);
-                                                            fhInvMassSumAll10MevBin3->Fill(fmassD0[0]);
-                                                            fhInvMassSumAll5MevBin3->Fill(fmassD0[0]);}
-        else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) {fhInvMassD0Sum10MevBin4->Fill(fmassD0[0]);
-                                                            fhInvMassD0Sum5MevBin4->Fill(fmassD0[0]);
-                                                            fhInvMassSumAll10MevBin4->Fill(fmassD0[0]);
-                                                            fhInvMassSumAll5MevBin4->Fill(fmassD0[0]);}
-        else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5])    {fhInvMassD0Sum10MevBin5->Fill(fmassD0[0]);
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) {fhInvMassD0Sum10MevBin2->Fill(fmassD0[0]);
+                                                             fhInvMassD0Sum5MevBin2->Fill(fmassD0[0]);
+                                                             fhInvMassSumAll10MevBin2->Fill(fmassD0[0]);
+                                                             fhInvMassSumAll5MevBin2->Fill(fmassD0[0]);}
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) {fhInvMassD0Sum10MevBin3->Fill(fmassD0[0]);
+                                                             fhInvMassD0Sum5MevBin3->Fill(fmassD0[0]);
+                                                             fhInvMassSumAll10MevBin3->Fill(fmassD0[0]);
+                                                             fhInvMassSumAll5MevBin3->Fill(fmassD0[0]);}
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) {fhInvMassD0Sum10MevBin4->Fill(fmassD0[0]);
+                                                             fhInvMassD0Sum5MevBin4->Fill(fmassD0[0]);
+                                                             fhInvMassSumAll10MevBin4->Fill(fmassD0[0]);
+                                                             fhInvMassSumAll5MevBin4->Fill(fmassD0[0]);}
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5])  {fhInvMassD0Sum10MevBin5->Fill(fmassD0[0]);
                                                                fhInvMassD0Sum5MevBin5->Fill(fmassD0[0]);
                                                                fhInvMassSumAll10MevBin5->Fill(fmassD0[0]);
                                                                fhInvMassSumAll5MevBin5->Fill(fmassD0[0]);} 
@@ -892,23 +930,23 @@ void AliAnalysisTaskSESelectHF4Prong::UserExec(Option_t */*option*/)
       if(hypD02==1) {
        fPtSel->Fill(ptPart);
         fScatterP4PID->Fill(trk3->Pt(),trk1->Pt());
-        if (ptPart < ptBinH[1])                         {fhInvMassD0Sum10MevBin1->Fill(fmassD0[1]);
+        if (ptPart < fPtBinH[1])                         {fhInvMassD0Sum10MevBin1->Fill(fmassD0[1]);
                                                          fhInvMassD0Sum5MevBin1->Fill(fmassD0[1]);
                                                          fhInvMassSumAll10MevBin1->Fill(fmassD0[1]);
                                                          fhInvMassSumAll5MevBin1->Fill(fmassD0[1]);}
-        else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) {fhInvMassD0Sum10MevBin2->Fill(fmassD0[1]);
-                                                            fhInvMassD0Sum5MevBin2->Fill(fmassD0[1]);
-                                                            fhInvMassSumAll10MevBin2->Fill(fmassD0[1]);
-                                                            fhInvMassSumAll5MevBin2->Fill(fmassD0[1]);}
-        else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) {fhInvMassD0Sum10MevBin3->Fill(fmassD0[1]);
-                                                            fhInvMassD0Sum5MevBin3->Fill(fmassD0[1]);
-                                                            fhInvMassSumAll10MevBin3->Fill(fmassD0[1]);
-                                                            fhInvMassSumAll5MevBin3->Fill(fmassD0[1]);}
-        else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) {fhInvMassD0Sum10MevBin4->Fill(fmassD0[1]);
-                                                            fhInvMassD0Sum5MevBin4->Fill(fmassD0[1]);
-                                                            fhInvMassSumAll10MevBin4->Fill(fmassD0[1]);
-                                                            fhInvMassSumAll5MevBin4->Fill(fmassD0[1]);}
-        else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5])    {fhInvMassD0Sum10MevBin5->Fill(fmassD0[1]);
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) {fhInvMassD0Sum10MevBin2->Fill(fmassD0[1]);
+                                                             fhInvMassD0Sum5MevBin2->Fill(fmassD0[1]);
+                                                             fhInvMassSumAll10MevBin2->Fill(fmassD0[1]);
+                                                             fhInvMassSumAll5MevBin2->Fill(fmassD0[1]);}
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) {fhInvMassD0Sum10MevBin3->Fill(fmassD0[1]);
+                                                             fhInvMassD0Sum5MevBin3->Fill(fmassD0[1]);
+                                                             fhInvMassSumAll10MevBin3->Fill(fmassD0[1]);
+                                                             fhInvMassSumAll5MevBin3->Fill(fmassD0[1]);}
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) {fhInvMassD0Sum10MevBin4->Fill(fmassD0[1]);
+                                                             fhInvMassD0Sum5MevBin4->Fill(fmassD0[1]);
+                                                             fhInvMassSumAll10MevBin4->Fill(fmassD0[1]);
+                                                             fhInvMassSumAll5MevBin4->Fill(fmassD0[1]);}
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5])  {fhInvMassD0Sum10MevBin5->Fill(fmassD0[1]);
                                                                 fhInvMassD0Sum5MevBin5->Fill(fmassD0[1]);
                                                                 fhInvMassSumAll10MevBin5->Fill(fmassD0[1]);
                                                                 fhInvMassSumAll5MevBin5->Fill(fmassD0[1]);}
@@ -917,23 +955,23 @@ void AliAnalysisTaskSESelectHF4Prong::UserExec(Option_t */*option*/)
       if(hypD0bar1==1) {
        fPtSel->Fill(ptPart);
         fScatterP4PID->Fill(trk0->Pt(),trk2->Pt());
-        if (ptPart < ptBinH[1])                         {fhInvMassD0barSum10MevBin1->Fill(fmassD0bar[0]);
+        if (ptPart < fPtBinH[1])                         {fhInvMassD0barSum10MevBin1->Fill(fmassD0bar[0]);
                                                          fhInvMassD0barSum5MevBin1->Fill(fmassD0bar[0]);
                                                          fhInvMassSumAll10MevBin1->Fill(fmassD0bar[0]);
                                                          fhInvMassSumAll5MevBin1->Fill(fmassD0bar[0]);}
-        else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) {fhInvMassD0barSum10MevBin2->Fill(fmassD0bar[0]);
-                                                            fhInvMassD0barSum5MevBin2->Fill(fmassD0bar[0]);
-                                                            fhInvMassSumAll10MevBin2->Fill(fmassD0bar[0]);
-                                                                    fhInvMassSumAll5MevBin2->Fill(fmassD0bar[0]);}
-        else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) {fhInvMassD0barSum10MevBin3->Fill(fmassD0bar[0]);
-                                                            fhInvMassD0barSum5MevBin3->Fill(fmassD0bar[0]);
-                                                            fhInvMassSumAll10MevBin3->Fill(fmassD0bar[0]);
-                                                            fhInvMassSumAll5MevBin3->Fill(fmassD0bar[0]);}
-        else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) {fhInvMassD0barSum10MevBin4->Fill(fmassD0bar[0]);
-                                                            fhInvMassD0barSum5MevBin4->Fill(fmassD0bar[0]);
-                                                            fhInvMassSumAll10MevBin4->Fill(fmassD0bar[0]);
-                                                            fhInvMassSumAll5MevBin4->Fill(fmassD0bar[0]);}
-        else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5])    {fhInvMassD0barSum10MevBin5->Fill(fmassD0bar[0]);
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) {fhInvMassD0barSum10MevBin2->Fill(fmassD0bar[0]);
+                                                             fhInvMassD0barSum5MevBin2->Fill(fmassD0bar[0]);
+                                                             fhInvMassSumAll10MevBin2->Fill(fmassD0bar[0]);
+                                                                     fhInvMassSumAll5MevBin2->Fill(fmassD0bar[0]);}
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) {fhInvMassD0barSum10MevBin3->Fill(fmassD0bar[0]);
+                                                             fhInvMassD0barSum5MevBin3->Fill(fmassD0bar[0]);
+                                                             fhInvMassSumAll10MevBin3->Fill(fmassD0bar[0]);
+                                                             fhInvMassSumAll5MevBin3->Fill(fmassD0bar[0]);}
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) {fhInvMassD0barSum10MevBin4->Fill(fmassD0bar[0]);
+                                                             fhInvMassD0barSum5MevBin4->Fill(fmassD0bar[0]);
+                                                             fhInvMassSumAll10MevBin4->Fill(fmassD0bar[0]);
+                                                             fhInvMassSumAll5MevBin4->Fill(fmassD0bar[0]);}
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5])  {fhInvMassD0barSum10MevBin5->Fill(fmassD0bar[0]);
                                                                fhInvMassD0barSum5MevBin5->Fill(fmassD0bar[0]);
                                                                fhInvMassSumAll10MevBin5->Fill(fmassD0bar[0]);
                                                                fhInvMassSumAll5MevBin5->Fill(fmassD0bar[0]);}   
@@ -942,29 +980,28 @@ void AliAnalysisTaskSESelectHF4Prong::UserExec(Option_t */*option*/)
       if(hypD0bar2==1) {
        fPtSel->Fill(ptPart);
         fScatterP4PID->Fill(trk2->Pt(),trk0->Pt());
-        if (ptPart < ptBinH[1])                         {fhInvMassD0barSum10MevBin1->Fill(fmassD0bar[1]);
+        if (ptPart < fPtBinH[1])                         {fhInvMassD0barSum10MevBin1->Fill(fmassD0bar[1]);
                                                          fhInvMassD0barSum5MevBin1->Fill(fmassD0bar[1]);
                                                          fhInvMassSumAll10MevBin1->Fill(fmassD0bar[1]);
                                                          fhInvMassSumAll5MevBin1->Fill(fmassD0bar[1]);}
-        else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) {fhInvMassD0barSum10MevBin2->Fill(fmassD0bar[1]);
-                                                            fhInvMassD0barSum5MevBin2->Fill(fmassD0bar[1]);
-                                                            fhInvMassSumAll10MevBin2->Fill(fmassD0bar[1]);
-                                                                    fhInvMassSumAll5MevBin2->Fill(fmassD0bar[1]);}
-        else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) {fhInvMassD0barSum10MevBin3->Fill(fmassD0bar[1]);
-                                                            fhInvMassD0barSum5MevBin3->Fill(fmassD0bar[1]);
-                                                            fhInvMassSumAll10MevBin3->Fill(fmassD0bar[1]);
-                                                            fhInvMassSumAll5MevBin3->Fill(fmassD0bar[1]);}
-        else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) {fhInvMassD0barSum10MevBin4->Fill(fmassD0bar[1]);
-                                                            fhInvMassD0barSum5MevBin4->Fill(fmassD0bar[1]);
-                                                            fhInvMassSumAll10MevBin4->Fill(fmassD0bar[1]);
-                                                            fhInvMassSumAll5MevBin4->Fill(fmassD0bar[1]);}
-        else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5])    {fhInvMassD0barSum10MevBin5->Fill(fmassD0bar[1]);
-                                                               fhInvMassD0barSum5MevBin5->Fill(fmassD0bar[1]);
-                                                               fhInvMassSumAll10MevBin5->Fill(fmassD0bar[1]);
-                                                               fhInvMassSumAll5MevBin5->Fill(fmassD0bar[1]);} 
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) {fhInvMassD0barSum10MevBin2->Fill(fmassD0bar[1]);
+                                                             fhInvMassD0barSum5MevBin2->Fill(fmassD0bar[1]);
+                                                             fhInvMassSumAll10MevBin2->Fill(fmassD0bar[1]);
+                                                                     fhInvMassSumAll5MevBin2->Fill(fmassD0bar[1]);}
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) {fhInvMassD0barSum10MevBin3->Fill(fmassD0bar[1]);
+                                                             fhInvMassD0barSum5MevBin3->Fill(fmassD0bar[1]);
+                                                             fhInvMassSumAll10MevBin3->Fill(fmassD0bar[1]);
+                                                             fhInvMassSumAll5MevBin3->Fill(fmassD0bar[1]);}
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) {fhInvMassD0barSum10MevBin4->Fill(fmassD0bar[1]);
+                                                             fhInvMassD0barSum5MevBin4->Fill(fmassD0bar[1]);
+                                                             fhInvMassSumAll10MevBin4->Fill(fmassD0bar[1]);
+                                                             fhInvMassSumAll5MevBin4->Fill(fmassD0bar[1]);}
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5])  {fhInvMassD0barSum10MevBin5->Fill(fmassD0bar[1]);
+                                                               fhInvMassD0barSum5MevBin5->Fill(fmassD0bar[1]);
+                                                               fhInvMassSumAll10MevBin5->Fill(fmassD0bar[1]);
+                                                               fhInvMassSumAll5MevBin5->Fill(fmassD0bar[1]);} 
      }
    }
-
       PostData(1,fOutput);
       PostData(2,fOutput2);
       PostData(3,fOutput3);
@@ -983,7 +1020,7 @@ void AliAnalysisTaskSESelectHF4Prong::UserExec(Option_t */*option*/)
       dOut->SetOwnPrimaryVtx((AliAODVertex*)((dIn->GetOwnPrimaryVtx())->Clone()));
       v->SetParent(dOut); 
       }
-    } //end of selection loop (starts with fSelected == 1, 2 or 3)
+    }
     if(unsetvtx) dIn->UnsetOwnPrimaryVtx();
   } // end loop on D0->K3pi
 
@@ -992,11 +1029,328 @@ void AliAnalysisTaskSESelectHF4Prong::UserExec(Option_t */*option*/)
   return;
 }
 
+//________________________________________________________________________
+void AliAnalysisTaskSESelectHF4Prong::AnalysisReflection(AliAODEvent* aodIn, AliAODRecoDecayHF4Prong* d, Int_t hypD01, Int_t hypD02, Int_t hypD0bar1, Int_t hypD0bar2)
+{
+  /*
+  ---STUDY OF REFLECTIONS ON CANDIDATE IN ANALYSIS---
+  Layers of TH2F fhRefl:
+  0 = all candidates, single hyps selected only;
+  1 = all candidates, multi hyps selected only;
+  2 = true D0/D0bar only, all hypotheses selected;
+  3 = true D0/D0bar only, single hyps selected only; -> I assume that the unique selected hypotesis is the right one!
+  4 = true D0/D0bar only, multi hyps selected only;
+  5 = true D0/D0bar only, multi TRUE hyps selected only;
+  6 = true D0/D0bar only, multi FAKE hyps selected only;
+  7 = false D0/D0bar only (background)
+  Layers of TH2F fhReflD0 (idem for D0bar, inverting particles):
+  0 = true D0 only, true hypothesis for both single and multi hyps cases;
+  1 = true D0 only, true hypothesis for single hyps case only; -> I assume that the unique selected hypotesis is the right one!
+  2 = true D0 only, true hypothesis for multi hyps case only;
+  3 = true D0 only, other D0 wrong hypothesis (multi case, obviously);
+  4 = true D0 only, D0bar1 wrong hypothesis (multi case, obviously);
+  5 = true D0 only, D0bar2 wrong hypothesis (multi case, obviously);
+  6 = true D0 only, D0bar1 + D0bar2 wrong hypotheses
+  */
+
+  Int_t flagMult = hypD01+hypD02+hypD0bar1+hypD0bar2; //single or multi hyps selected
+  Int_t flagLayer1 = -1, flagLayer2 = -1, flagLayer3 = -1, flagLayer4 = -1; //to select layers in fhRefl
+
+  if (flagMult == 1) {flagLayer1 = 0; flagLayer2 = 0; flagLayer3 = 0; flagLayer4 = 0;}
+  else if (flagMult == 2 || flagMult == 3 || flagMult == 4) {flagLayer1 = 1; flagLayer2 = 1; flagLayer3 = 1; flagLayer4 = 1;}
+
+  FillReflHistos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, flagLayer1, flagLayer2, flagLayer3, flagLayer4);  //Fill layers 0 and 1
+
+  if (fMCTruth==0) return;
+
+  //start of MC Truth phase
+  TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodIn->FindListObject(AliAODMCParticle::StdBranchName()));
+  if (!mcArray) {
+    AliError("Could not find Monte-Carlo in AOD");
+    return;}
+
+  Int_t pdgCand = 421;
+  Int_t pdgDgCharm4Prong[4]={321,211,211,211}; //pdg of daughters
+
+  Int_t mcLabel = d->MatchToMC(pdgCand,mcArray,4,pdgDgCharm4Prong); //selection of true or false candidate (regardless of hypothesis) through MCtruth
+  printf("MatchToMC = %d\n",mcLabel); 
+
+  if (mcLabel==-1) { //fill layer 7 (background)
+  FillReflHistos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, 7, 7, 7, 7);
+  return;}
+
+  FillReflHistos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, 2, 2, 2, 2); //fill layer 2 (true D0/D0bar)
+
+  Int_t truthHyp = StudyMCTruth(mcArray, d); //calls function which studies which hypothesis is true for candidate
+
+  if (flagMult == 1) { //fill layer 3 (true D0/D0bar - single hyps only)
+    FillReflHistos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, 3, 3, 3, 3);
+    switch(truthHyp) { //fills fhReflD0 and fhReflD0bar with layers containing all true D0 or D0bars, for single and all hypotheses (layers 1, 0)
+      case(1): FillReflD0Histos(d, hypD01, 0, 0, 0, 0, 0, 0, 0);
+               FillReflD0Histos(d, hypD01, 0, 0, 0, 1, 1, 1, 1); //here I discard the cases in which only a wrong hyp is selected (very very few)
+               break;
+      case(2): FillReflD0Histos(d, 0, hypD02, 0, 0, 0, 0, 0, 0);
+               FillReflD0Histos(d, 0, hypD02, 0, 0, 1, 1, 1, 1);
+               break;
+      case(3): FillReflD0barHistos(d, 0, 0, hypD0bar1, 0, 0, 0, 0, 0); 
+              FillReflD0barHistos(d, 0, 0, hypD0bar1, 0, 1, 1, 1, 1);
+              break;
+      case(4): FillReflD0barHistos(d, 0, 0, 0, hypD0bar2, 0, 0, 0, 0); 
+              FillReflD0barHistos(d, 0, 0, 0, hypD0bar2, 1, 1, 1, 1);
+              break;
+    }
+  }
+  else { 
+    flagLayer1 = 6; flagLayer2 = 6; flagLayer3 = 6; flagLayer4 = 6; 
+    switch(truthHyp) { //fills fhReflD0 and fhReflD0bar with layers containing all true D0 or D0bars, for multi and all hypotheses (layers 2, 0)
+      case(1): FillReflD0Histos(d, hypD01, 0, 0, 0, 0, 0, 0, 0);
+               FillReflD0Histos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, 2, 3, 4, 5);
+               FillReflD0Histos(d, 0, 0, hypD0bar1, hypD0bar2, 2, 3, 6, 6); //merge of opposite particle hyps (D0bar) in layer 6
+              flagLayer1 = 5;
+               break;
+      case(2): FillReflD0Histos(d, 0, hypD02, 0, 0, 0, 0, 0, 0);
+               FillReflD0Histos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, 3, 2, 4, 5);
+               FillReflD0Histos(d, 0, 0, hypD0bar1, hypD0bar2, 0, 0, 6, 6); //merge of opposite particle hyps (D0bar) in layer 6
+              flagLayer2 = 5;
+               break;
+      case(3): FillReflD0barHistos(d, 0, 0, hypD0bar1, 0, 0, 0, 0, 0); 
+              FillReflD0barHistos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, 4, 5, 2, 3);
+               FillReflD0barHistos(d, hypD01, hypD02, 0, 0, 6, 6, 0, 0); //merge of opposite particle hyps (D0) in layer 6
+              flagLayer3 = 5;
+              break;
+      case(4): FillReflD0barHistos(d, 0, 0, 0, hypD0bar2, 0, 0, 0, 0); 
+              FillReflD0barHistos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, 4, 5, 3, 2);
+               FillReflD0barHistos(d, hypD01, hypD02, 0, 0, 6, 6, 0, 0); //merge of opposite particle hyps (D0) in layer 6
+              flagLayer4 = 5;
+              break;
+    }
+    FillReflHistos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, flagLayer1, flagLayer2, flagLayer3, flagLayer4); //fill layers 5 and 6 (true and false hyps for multi)
+    FillReflHistos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, 4, 4, 4, 4); //fill layer 4 (true D0/D0bar - multi hyps only)
+  }
+}
+
+//________________________________________________________________________
+Int_t AliAnalysisTaskSESelectHF4Prong::StudyMCTruth(TClonesArray* mcArray, AliAODRecoDecayHF4Prong* d)
+{
+  /* 
+  ---STUDY OF MCTRUTH ON CANDIDATE IN ANALYSIS---
+  Flag Truth (output):
+  0 = problems in daughter tracks found
+  1 = candidate is D01 (piKpipi)
+  2 = candidate is D02 (pipipiK)
+  3 = candidate is D0bar1 (Kpipipi)
+  4 = candidate is D0bar2 (pipiKpi)
+  */
+
+  Int_t truthHyp = 0;
+
+  AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
+  AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1);
+  AliAODTrack *trk2 = (AliAODTrack*)d->GetDaughter(2);
+  AliAODTrack *trk3 = (AliAODTrack*)d->GetDaughter(3);
+  Int_t labels[4];
+  Int_t pdg[4];
+  labels[0] = trk0->GetLabel();
+  labels[1] = trk1->GetLabel();
+  labels[2] = trk2->GetLabel();
+  labels[3] = trk3->GetLabel();
+  if (labels[0]<=0 || labels[1]<=0 || labels[2]<=0 || labels[3]<=0) {AliWarning("Negative Label for daughter, skipping"); return truthHyp;}
+  AliAODMCParticle* mc0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labels[0]));
+  AliAODMCParticle* mc1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labels[1]));
+  AliAODMCParticle* mc2 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labels[2]));
+  AliAODMCParticle* mc3 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labels[3]));
+  if (!mc0 || !mc1 || !mc2 || !mc3) {AliWarning("At least one Daughter Particle not found in tree, skipping"); return truthHyp;}
+  pdg[0] = TMath::Abs(mc0->GetPdgCode());
+  pdg[1] = TMath::Abs(mc1->GetPdgCode());
+  pdg[2] = TMath::Abs(mc2->GetPdgCode());
+  pdg[3] = TMath::Abs(mc3->GetPdgCode());
+  if (pdg[0]==211 && pdg[1]==321 && pdg[2]==211 && pdg[3]==211) truthHyp = 1;
+  if (pdg[0]==211 && pdg[1]==211 && pdg[2]==211 && pdg[3]==321) truthHyp = 2;
+  if (pdg[0]==321 && pdg[1]==211 && pdg[2]==211 && pdg[3]==211) truthHyp = 3;
+  if (pdg[0]==211 && pdg[1]==211 && pdg[2]==321 && pdg[3]==211) truthHyp = 4;
+
+  return truthHyp;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSESelectHF4Prong::FillReflHistos(AliAODRecoDecayHF4Prong* d, Int_t hypD01, Int_t hypD02, Int_t hypD0bar1, Int_t hypD0bar2, Int_t flagLayer1, Int_t flagLayer2, Int_t flagLayer3, Int_t flagLayer4)
+{
+
+  Double_t ptPart = d->Pt();
+
+    if (ptPart > fPtBinH[0]) {
+      // D01 hyp.
+      if(hypD01==1) {
+        if (ptPart < fPtBinH[1])                             fhReflBin1->Fill(fmassD0[0],(double)flagLayer1);
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflBin2->Fill(fmassD0[0],(double)flagLayer1);
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflBin3->Fill(fmassD0[0],(double)flagLayer1);
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflBin4->Fill(fmassD0[0],(double)flagLayer1);
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflBin5->Fill(fmassD0[0],(double)flagLayer1);
+      }
+      // D02 hyp.
+      if(hypD02==1) {
+        if (ptPart < fPtBinH[1])                             fhReflBin1->Fill(fmassD0[1],(double)flagLayer2);
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflBin2->Fill(fmassD0[1],(double)flagLayer2);
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflBin3->Fill(fmassD0[1],(double)flagLayer2);
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflBin4->Fill(fmassD0[1],(double)flagLayer2);
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflBin5->Fill(fmassD0[1],(double)flagLayer2);
+      }
+      // D0bar1 hyp.
+      if(hypD0bar1==1) {
+       if (ptPart < fPtBinH[1])                              fhReflBin1->Fill(fmassD0bar[0],(double)flagLayer3);
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflBin2->Fill(fmassD0bar[0],(double)flagLayer3);
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflBin3->Fill(fmassD0bar[0],(double)flagLayer3);
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflBin4->Fill(fmassD0bar[0],(double)flagLayer3);
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflBin5->Fill(fmassD0bar[0],(double)flagLayer3);
+      } 
+      // D0bar2 hyp.
+      if(hypD0bar2==1) {
+        if (ptPart < fPtBinH[1])                             fhReflBin1->Fill(fmassD0bar[1],(double)flagLayer4);
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflBin2->Fill(fmassD0bar[1],(double)flagLayer4);
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflBin3->Fill(fmassD0bar[1],(double)flagLayer4);
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflBin4->Fill(fmassD0bar[1],(double)flagLayer4);
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflBin5->Fill(fmassD0bar[1],(double)flagLayer4);
+      }
+    }
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSESelectHF4Prong::FillReflD0Histos(AliAODRecoDecayHF4Prong* d, Int_t hypD01, Int_t hypD02, Int_t hypD0bar1, Int_t hypD0bar2, Int_t flagLayforD01, Int_t flagLayforD02, Int_t flagLayforD03, Int_t flagLayforD04)
+{
+
+  Double_t ptPart = d->Pt();
+
+    if (ptPart > fPtBinH[0]) {
+      // D01 hyp.
+      if(hypD01==1) {
+        if (ptPart < fPtBinH[1])                             fhReflD0Bin1->Fill(fmassD0[0],(double)flagLayforD01);
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflD0Bin2->Fill(fmassD0[0],(double)flagLayforD01);
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflD0Bin3->Fill(fmassD0[0],(double)flagLayforD01);
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflD0Bin4->Fill(fmassD0[0],(double)flagLayforD01);
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflD0Bin5->Fill(fmassD0[0],(double)flagLayforD01);
+      }
+      // D02 hyp.
+      if(hypD02==1) {
+        if (ptPart < fPtBinH[1])                             fhReflD0Bin1->Fill(fmassD0[1],(double)flagLayforD02);
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflD0Bin2->Fill(fmassD0[1],(double)flagLayforD02);
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflD0Bin3->Fill(fmassD0[1],(double)flagLayforD02);
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflD0Bin4->Fill(fmassD0[1],(double)flagLayforD02);
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflD0Bin5->Fill(fmassD0[1],(double)flagLayforD02);
+      }
+      // D0bar1 hyp.
+      if(hypD0bar1==1) {
+       if (ptPart < fPtBinH[1])                              fhReflD0Bin1->Fill(fmassD0bar[0],(double)flagLayforD03);
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflD0Bin2->Fill(fmassD0bar[0],(double)flagLayforD03);
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflD0Bin3->Fill(fmassD0bar[0],(double)flagLayforD03);
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflD0Bin4->Fill(fmassD0bar[0],(double)flagLayforD03);
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflD0Bin5->Fill(fmassD0bar[0],(double)flagLayforD03);
+      } 
+      // D0bar2 hyp.
+      if(hypD0bar2==1) {
+        if (ptPart < fPtBinH[1])                             fhReflD0Bin1->Fill(fmassD0bar[1],(double)flagLayforD04);
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflD0Bin2->Fill(fmassD0bar[1],(double)flagLayforD04);
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflD0Bin3->Fill(fmassD0bar[1],(double)flagLayforD04);
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflD0Bin4->Fill(fmassD0bar[1],(double)flagLayforD04);
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflD0Bin5->Fill(fmassD0bar[1],(double)flagLayforD04);
+      }
+    }
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSESelectHF4Prong::FillReflD0barHistos(AliAODRecoDecayHF4Prong* d, Int_t hypD01, Int_t hypD02, Int_t hypD0bar1, Int_t hypD0bar2, Int_t flagLayforD0bar1, Int_t flagLayforD0bar2, Int_t flagLayforD0bar3, Int_t flagLayforD0bar4)
+{
+
+  Double_t ptPart = d->Pt();
+
+    if (ptPart > fPtBinH[0]) {
+      // D01 hyp.
+      if(hypD01==1) {
+        if (ptPart < fPtBinH[1])                             fhReflD0barBin1->Fill(fmassD0[0],(double)flagLayforD0bar1);
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflD0barBin2->Fill(fmassD0[0],(double)flagLayforD0bar1);
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflD0barBin3->Fill(fmassD0[0],(double)flagLayforD0bar1);
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflD0barBin4->Fill(fmassD0[0],(double)flagLayforD0bar1);
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflD0barBin5->Fill(fmassD0[0],(double)flagLayforD0bar1);
+      }
+      // D02 hyp.
+      if(hypD02==1) {
+        if (ptPart < fPtBinH[1])                             fhReflD0barBin1->Fill(fmassD0[1],(double)flagLayforD0bar2);
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflD0barBin2->Fill(fmassD0[1],(double)flagLayforD0bar2);
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflD0barBin3->Fill(fmassD0[1],(double)flagLayforD0bar2);
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflD0barBin4->Fill(fmassD0[1],(double)flagLayforD0bar2);
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflD0barBin5->Fill(fmassD0[1],(double)flagLayforD0bar2);
+      }
+      // D0bar1 hyp.
+      if(hypD0bar1==1) {
+       if (ptPart < fPtBinH[1])                              fhReflD0barBin1->Fill(fmassD0bar[0],(double)flagLayforD0bar3);
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflD0barBin2->Fill(fmassD0bar[0],(double)flagLayforD0bar3);
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflD0barBin3->Fill(fmassD0bar[0],(double)flagLayforD0bar3);
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflD0barBin4->Fill(fmassD0bar[0],(double)flagLayforD0bar3);
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflD0barBin5->Fill(fmassD0bar[0],(double)flagLayforD0bar3);
+      } 
+      // D0bar2 hyp.
+      if(hypD0bar2==1) {
+        if (ptPart < fPtBinH[1])                             fhReflD0barBin1->Fill(fmassD0bar[1],(double)flagLayforD0bar4);
+        else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflD0barBin2->Fill(fmassD0bar[1],(double)flagLayforD0bar4);
+        else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflD0barBin3->Fill(fmassD0bar[1],(double)flagLayforD0bar4);
+        else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflD0barBin4->Fill(fmassD0bar[1],(double)flagLayforD0bar4);
+        else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflD0barBin5->Fill(fmassD0bar[1],(double)flagLayforD0bar4);
+      }
+    }
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSESelectHF4Prong::SetPtBinH(Double_t* ptlimits)
+{
+  for(int i=0; i<6; i++) {fPtBinH[i]=ptlimits[i];}
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSESelectHF4Prong::PrintPtBinHandMCFlag()
+{
+  printf("PtBin limits---------\n");
+  for (int i=0; i<5; i++) {
+    printf("Bin %d = %.1f to %.1f\n",i+1,fPtBinH[i],fPtBinH[i+1]);
+  }
+  printf("MC Truth = %d\n",fMCTruth);
+  printf("---------------------\n");
+}
+
 //________________________________________________________________________
 void AliAnalysisTaskSESelectHF4Prong::Terminate(Option_t */*option*/)
 {
   // Terminate analysis
   //
+
+/*
+  Double_t entries[10] = {0,0,0,0,0,0,0,0,0,0};
+  for(int i=1;i<=fhReflD0Bin1->GetNbinsX();i++) {
+    for(int j=1;j<=fhReflD0Bin1->GetNbinsY();j++) {
+      entries[0] += fhReflD0Bin1->GetBinContent(i,j);
+      entries[1] += fhReflD0Bin2->GetBinContent(i,j);
+      entries[2] += fhReflD0Bin3->GetBinContent(i,j);
+      entries[3] += fhReflD0Bin4->GetBinContent(i,j);
+      entries[4] += fhReflD0Bin5->GetBinContent(i,j);
+      entries[5] += fhReflD0barBin1->GetBinContent(i,j);
+      entries[6] += fhReflD0barBin2->GetBinContent(i,j);
+      entries[7] += fhReflD0barBin3->GetBinContent(i,j);
+      entries[8] += fhReflD0barBin4->GetBinContent(i,j);
+      entries[9] += fhReflD0barBin5->GetBinContent(i,j);
+    }
+  }
+  fhReflD0Bin1->SetEntries(entries[0]);
+  fhReflD0Bin2->SetEntries(entries[1]);  
+  fhReflD0Bin3->SetEntries(entries[2]);
+  fhReflD0Bin4->SetEntries(entries[3]);
+  fhReflD0Bin5->SetEntries(entries[4]);
+  fhReflD0barBin1->SetEntries(entries[5]);
+  fhReflD0barBin2->SetEntries(entries[6]);
+  fhReflD0barBin3->SetEntries(entries[7]);
+  fhReflD0barBin4->SetEntries(entries[8]);
+  fhReflD0barBin5->SetEntries(entries[9]);*/
+
   if(fDebug > 1) printf("AnalysisTaskSESelectHF4Prong: Terminate() \n");
 
   fOutput = dynamic_cast<TList*> (GetOutputData(1));
@@ -1019,11 +1373,12 @@ void AliAnalysisTaskSESelectHF4Prong::Terminate(Option_t */*option*/)
     printf("ERROR: fOutput not available\n");
     return;
   }
+  fOutput5 = dynamic_cast<TList*> (GetOutputData(5));
   if (!fOutput5) {
     printf("ERROR: fOutput not available\n");
     return;
   }
-  fOutputC = dynamic_cast<TList*> (GetOutputData(5));
+  fOutputC = dynamic_cast<TList*> (GetOutputData(6));
   if (!fOutputC) {
     printf("ERROR: fOutputC not available\n");
     return;