--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/////////////////////////////////////////////////////////////
+//
+// AliAnalysisTaskSE for the selection of heavy flavor
+// decay candidates and creation a stand-alone AOD for
+// 4prong D0 decay.
+//
+// Author: A.Dainese, andrea.dainese@lnl.infn.it
+// F.Colamaria, fabio.colamaria@ba.infn.it
+/////////////////////////////////////////////////////////////
+
+
+#include "Riostream.h"
+#include "TFile.h"
+#include "TList.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TClonesArray.h"
+#include "TDatabasePDG.h"
+#include "TROOT.h"
+#include "TCanvas.h"
+#include "TBits.h"
+#include "TNtuple.h"
+
+#include "AliAnalysisDataSlot.h"
+#include "AliAnalysisDataContainer.h"
+#include "AliAnalysisManager.h"
+#include "AliAODHandler.h"
+#include "AliAODEvent.h"
+#include "AliAODVertex.h"
+#include "AliAODTrack.h"
+#include "AliAODRecoDecayHF4Prong.h"
+#include "AliAnalysisVertexingHF.h"
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisTaskSESelectHF4Prong.h"
+#include "AliAODPidHF.h"
+#include "AliRDHFCuts.h"
+
+ClassImp(AliAnalysisTaskSESelectHF4Prong)
+
+
+//________________________________________________________________________
+AliAnalysisTaskSESelectHF4Prong::AliAnalysisTaskSESelectHF4Prong():
+AliAnalysisTaskSE(),
+fVerticesHFTClArr(0),
+fCharm4ProngTClArr(0),
+fOutput(0),
+fOutput2(0),
+fOutput3(0),
+fOutput4(0),
+fOutput5(0),
+fOutputC(0),
+fhInvMassD0Sum_10Mev_Bin1(0),
+fhInvMassD0barSum_10Mev_Bin1(0),
+fhInvMassSumAll_10Mev_Bin1(0),
+fhInvMassD0Sum_5Mev_Bin1(0),
+fhInvMassD0barSum_5Mev_Bin1(0),
+fhInvMassSumAll_5Mev_Bin1(0),
+fhInvMassD0Sum_10Mev_Bin2(0),
+fhInvMassD0barSum_10Mev_Bin2(0),
+fhInvMassSumAll_10Mev_Bin2(0),
+fhInvMassD0Sum_5Mev_Bin2(0),
+fhInvMassD0barSum_5Mev_Bin2(0),
+fhInvMassSumAll_5Mev_Bin2(0),
+fhInvMassD0Sum_10Mev_Bin3(0),
+fhInvMassD0barSum_10Mev_Bin3(0),
+fhInvMassSumAll_10Mev_Bin3(0),
+fhInvMassD0Sum_5Mev_Bin3(0),
+fhInvMassD0barSum_5Mev_Bin3(0),
+fhInvMassSumAll_5Mev_Bin3(0),
+fhInvMassD0Sum_10Mev_Bin4(0),
+fhInvMassD0barSum_10Mev_Bin4(0),
+fhInvMassSumAll_10Mev_Bin4(0),
+fhInvMassD0Sum_5Mev_Bin4(0),
+fhInvMassD0barSum_5Mev_Bin4(0),
+fhInvMassSumAll_5Mev_Bin4(0),
+fhInvMassD0Sum_10Mev_Bin5(0),
+fhInvMassD0barSum_10Mev_Bin5(0),
+fhInvMassSumAll_10Mev_Bin5(0),
+fhInvMassD0Sum_5Mev_Bin5(0),
+fhInvMassD0barSum_5Mev_Bin5(0),
+fhInvMassSumAll_5Mev_Bin5(0),
+fhInvMassMultipleOnly_Bin1(0),
+fhInvMassMultipleOnly_Bin2(0),
+fhInvMassMultipleOnly_Bin3(0),
+fhInvMassMultipleOnly_Bin4(0),
+fhInvMassMultipleOnly_Bin5(0),
+fScatterP4PID(0),
+fPtVsY(0),
+fPtVsYAll(0),
+fEventCounter(0),
+fCutDCA(0),
+fCutDCA3(0),
+fCutDCA2(0),
+fCutDCA5(0),
+fCutVertexDist2(0),
+fCutVertexDist3(0),
+fCutVertexDist4(0),
+fCutCosinePoint(0),
+fCutPt(0),
+fCutY(0),
+fPIDSel(0),
+fPIDSel_Bin1(0),
+fPIDSel_Bin2(0),
+fPIDSel_Bin3(0),
+fPIDSel_Bin4(0),
+fPIDSel_Bin5(0),
+fMultipleHyps(0),
+fMultipleHypsType(0),
+fPtSel(0),
+fCuts(0),
+fVHF(0)
+{
+ // Default constructor
+}
+
+//________________________________________________________________________
+AliAnalysisTaskSESelectHF4Prong::AliAnalysisTaskSESelectHF4Prong(const char *name,AliRDHFCutsD0toKpipipi* cuts):
+AliAnalysisTaskSE(name),
+fVerticesHFTClArr(0),
+fCharm4ProngTClArr(0),
+fOutput(0),
+fOutput2(0),
+fOutput3(0),
+fOutput4(0),
+fOutput5(0),
+fOutputC(0),
+fhInvMassD0Sum_10Mev_Bin1(0),
+fhInvMassD0barSum_10Mev_Bin1(0),
+fhInvMassSumAll_10Mev_Bin1(0),
+fhInvMassD0Sum_5Mev_Bin1(0),
+fhInvMassD0barSum_5Mev_Bin1(0),
+fhInvMassSumAll_5Mev_Bin1(0),
+fhInvMassD0Sum_10Mev_Bin2(0),
+fhInvMassD0barSum_10Mev_Bin2(0),
+fhInvMassSumAll_10Mev_Bin2(0),
+fhInvMassD0Sum_5Mev_Bin2(0),
+fhInvMassD0barSum_5Mev_Bin2(0),
+fhInvMassSumAll_5Mev_Bin2(0),
+fhInvMassD0Sum_10Mev_Bin3(0),
+fhInvMassD0barSum_10Mev_Bin3(0),
+fhInvMassSumAll_10Mev_Bin3(0),
+fhInvMassD0Sum_5Mev_Bin3(0),
+fhInvMassD0barSum_5Mev_Bin3(0),
+fhInvMassSumAll_5Mev_Bin3(0),
+fhInvMassD0Sum_10Mev_Bin4(0),
+fhInvMassD0barSum_10Mev_Bin4(0),
+fhInvMassSumAll_10Mev_Bin4(0),
+fhInvMassD0Sum_5Mev_Bin4(0),
+fhInvMassD0barSum_5Mev_Bin4(0),
+fhInvMassSumAll_5Mev_Bin4(0),
+fhInvMassD0Sum_10Mev_Bin5(0),
+fhInvMassD0barSum_10Mev_Bin5(0),
+fhInvMassSumAll_10Mev_Bin5(0),
+fhInvMassD0Sum_5Mev_Bin5(0),
+fhInvMassD0barSum_5Mev_Bin5(0),
+fhInvMassSumAll_5Mev_Bin5(0),
+fhInvMassMultipleOnly_Bin1(0),
+fhInvMassMultipleOnly_Bin2(0),
+fhInvMassMultipleOnly_Bin3(0),
+fhInvMassMultipleOnly_Bin4(0),
+fhInvMassMultipleOnly_Bin5(0),
+fScatterP4PID(0),
+fPtVsY(0),
+fPtVsYAll(0),
+fEventCounter(0),
+fCutDCA(0),
+fCutDCA3(0),
+fCutDCA2(0),
+fCutDCA5(0),
+fCutVertexDist2(0),
+fCutVertexDist3(0),
+fCutVertexDist4(0),
+fCutCosinePoint(0),
+fCutPt(0),
+fCutY(0),
+fPIDSel(0),
+fPIDSel_Bin1(0),
+fPIDSel_Bin2(0),
+fPIDSel_Bin3(0),
+fPIDSel_Bin4(0),
+fPIDSel_Bin5(0),
+fMultipleHyps(0),
+fMultipleHypsType(0),
+fPtSel(0),
+fCuts(0),
+fVHF(0)
+{
+ // Standard constructor
+
+ fCuts=cuts;
+
+ // Input slot #0 works with an Ntuple
+// DefineInput(0, TTree::Class());
+
+ // Output slot #0 writes into a TTree container
+ // Output slots #1-6 writes into a TList container
+ DefineOutput(0, TTree::Class()); //default
+ DefineOutput(1, TList::Class()); //histos inv. mass bin1
+ DefineOutput(2, TList::Class()); //histos inv. mass bin2
+ DefineOutput(3, TList::Class()); //histos inv. mass bin3
+ DefineOutput(4, TList::Class()); //histos inv. mass bin4
+ DefineOutput(5, TList::Class()); //histos inv. mass bin5
+ DefineOutput(6, TList::Class()); //histos of cuts
+ DefineOutput(7, AliRDHFCutsD0toKpipipi::Class()); //cuts
+}
+
+//________________________________________________________________________
+AliAnalysisTaskSESelectHF4Prong::~AliAnalysisTaskSESelectHF4Prong()
+{
+ // Destructor
+
+ if (fOutput) {
+ delete fOutput;
+ fOutput = 0;
+ }
+ if (fOutput2) {
+ delete fOutput2;
+ fOutput2 = 0;
+ }
+ if (fOutput3) {
+ delete fOutput3;
+ fOutput3 = 0;
+ }
+ if (fOutput4) {
+ delete fOutput4;
+ fOutput4 = 0;
+ }
+ if (fOutput5) {
+ delete fOutput5;
+ fOutput5 = 0;
+ }
+ if (fOutputC) {
+ delete fOutputC;
+ fOutputC = 0;
+ }
+ if (fCuts) {
+ delete fCuts;
+ fCuts = 0;
+ }
+ if (fVHF) {
+ delete fVHF;
+ fVHF = 0;
+ }
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSESelectHF4Prong::Init()
+{
+ // Initialization
+
+ if(fDebug > 1) printf("AnalysisTaskSESelectHF4Prong::Init() \n");
+
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSESelectHF4Prong::UserCreateOutputObjects()
+{
+ // Create the output container
+ //
+ if(fDebug > 1) printf("AnalysisTaskSESelectHF4Prong::UserCreateOutputObjects() \n");
+
+ fVerticesHFTClArr = new TClonesArray("AliAODVertex", 0);
+ fVerticesHFTClArr->SetName("VerticesHF");
+ AddAODBranch("TClonesArray", &fVerticesHFTClArr);
+
+ fCharm4ProngTClArr = new TClonesArray("AliAODRecoDecayHF4Prong", 0);
+ fCharm4ProngTClArr->SetName("Charm4Prong");
+ AddAODBranch("TClonesArray", &fCharm4ProngTClArr);
+
+ fOutput = new TList();
+ fOutput->SetOwner();
+
+ fOutput2 = new TList();
+ fOutput2->SetOwner();
+
+ fOutput3 = new TList();
+ fOutput3->SetOwner();
+
+ fOutput4 = new TList();
+ fOutput4->SetOwner();
+
+ fOutput5 = new TList();
+ fOutput5->SetOwner();
+
+ fOutputC = new TList();
+ fOutputC->SetOwner();
+
+ fhInvMassD0Sum_10Mev_Bin1 = new TH1F("fhInvMassD0Sum_10Mev_Bin1", "D0 invariant mass Bin1 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
+ fhInvMassD0Sum_10Mev_Bin1->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0Sum_10Mev_Bin1->SetMinimum(0);
+ fOutput->Add(fhInvMassD0Sum_10Mev_Bin1);
+
+ fhInvMassD0barSum_10Mev_Bin1 = new TH1F("fhInvMassD0barSum_10Mev_Bin1", "D0bar invariant mass Bin1 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
+ fhInvMassD0barSum_10Mev_Bin1->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0barSum_10Mev_Bin1->SetMinimum(0);
+ fOutput->Add(fhInvMassD0barSum_10Mev_Bin1);
+
+ fhInvMassSumAll_10Mev_Bin1 = new TH1F("fhInvMassSumAll_10Mev_Bin1", "D0/D0bar invariant mass Bin1 (good hyps); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
+ fhInvMassSumAll_10Mev_Bin1->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassSumAll_10Mev_Bin1->SetMinimum(0);
+ fOutput->Add(fhInvMassSumAll_10Mev_Bin1);
+
+ fhInvMassD0Sum_5Mev_Bin1 = new TH1F("fhInvMassD0Sum_5Mev_Bin1", "D0 invariant mass Bin1 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
+ fhInvMassD0Sum_5Mev_Bin1->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0Sum_5Mev_Bin1->SetMinimum(0);
+ fOutput->Add(fhInvMassD0Sum_5Mev_Bin1);
+
+ fhInvMassD0barSum_5Mev_Bin1 = new TH1F("fhInvMassD0barSum_5Mev_Bin1", "D0bar invariant mass Bin1 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
+ fhInvMassD0barSum_5Mev_Bin1->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0barSum_5Mev_Bin1->SetMinimum(0);
+ fOutput->Add(fhInvMassD0barSum_5Mev_Bin1);
+
+ fhInvMassSumAll_5Mev_Bin1 = new TH1F("fhInvMassSumAll_5Mev_Bin1", "D0/D0bar invariant mass Bin1 (good hyps); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
+ fhInvMassSumAll_5Mev_Bin1->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassSumAll_5Mev_Bin1->SetMinimum(0);
+ fOutput->Add(fhInvMassSumAll_5Mev_Bin1);
+
+ fhInvMassD0Sum_10Mev_Bin2 = new TH1F("fhInvMassD0Sum_10Mev_Bin2", "D0 invariant mass Bin2 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
+ fhInvMassD0Sum_10Mev_Bin2->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0Sum_10Mev_Bin2->SetMinimum(0);
+ fOutput2->Add(fhInvMassD0Sum_10Mev_Bin2);
+
+ fhInvMassD0barSum_10Mev_Bin2 = new TH1F("fhInvMassD0barSum_10Mev_Bin2", "D0bar invariant mass Bin2 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
+ fhInvMassD0barSum_10Mev_Bin2->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0barSum_10Mev_Bin2->SetMinimum(0);
+ fOutput2->Add(fhInvMassD0barSum_10Mev_Bin2);
+
+ fhInvMassSumAll_10Mev_Bin2 = new TH1F("fhInvMassSumAll_10Mev_Bin2", "D0/D0bar invariant mass Bin2 (good hyps); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
+ fhInvMassSumAll_10Mev_Bin2->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassSumAll_10Mev_Bin2->SetMinimum(0);
+ fOutput2->Add(fhInvMassSumAll_10Mev_Bin2);
+
+ fhInvMassD0Sum_5Mev_Bin2 = new TH1F("fhInvMassD0Sum_5Mev_Bin2", "D0 invariant mass Bin2 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
+ fhInvMassD0Sum_5Mev_Bin2->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0Sum_5Mev_Bin2->SetMinimum(0);
+ fOutput2->Add(fhInvMassD0Sum_5Mev_Bin2);
+
+ fhInvMassD0barSum_5Mev_Bin2 = new TH1F("fhInvMassD0barSum_5Mev_Bin2", "D0bar invariant mass Bin2 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
+ fhInvMassD0barSum_5Mev_Bin2->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0barSum_5Mev_Bin2->SetMinimum(0);
+ fOutput2->Add(fhInvMassD0barSum_5Mev_Bin2);
+
+ fhInvMassSumAll_5Mev_Bin2 = new TH1F("fhInvMassSumAll_5Mev_Bin2", "D0/D0bar invariant mass Bin2 (good hyps); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
+ fhInvMassSumAll_5Mev_Bin2->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassSumAll_5Mev_Bin2->SetMinimum(0);
+ fOutput2->Add(fhInvMassSumAll_5Mev_Bin2);
+
+ fhInvMassD0Sum_10Mev_Bin3 = new TH1F("fhInvMassD0Sum_10Mev_Bin3", "D0 invariant mass Bin3 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
+ fhInvMassD0Sum_10Mev_Bin3->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0Sum_10Mev_Bin3->SetMinimum(0);
+ fOutput3->Add(fhInvMassD0Sum_10Mev_Bin3);
+
+ fhInvMassD0barSum_10Mev_Bin3 = new TH1F("fhInvMassD0barSum_10Mev_Bin3", "D0bar invariant mass Bin3 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
+ fhInvMassD0barSum_10Mev_Bin3->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0barSum_10Mev_Bin3->SetMinimum(0);
+ fOutput3->Add(fhInvMassD0barSum_10Mev_Bin3);
+
+ fhInvMassSumAll_10Mev_Bin3 = new TH1F("fhInvMassSumAll_10Mev_Bin3", "D0/D0bar invariant mass Bin3 (good hyps); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
+ fhInvMassSumAll_10Mev_Bin3->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassSumAll_10Mev_Bin3->SetMinimum(0);
+ fOutput3->Add(fhInvMassSumAll_10Mev_Bin3);
+
+ fhInvMassD0Sum_5Mev_Bin3 = new TH1F("fhInvMassD0Sum_5Mev_Bin3", "D0 invariant mass Bin3 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
+ fhInvMassD0Sum_5Mev_Bin3->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0Sum_5Mev_Bin3->SetMinimum(0);
+ fOutput3->Add(fhInvMassD0Sum_5Mev_Bin3);
+
+ fhInvMassD0barSum_5Mev_Bin3 = new TH1F("fhInvMassD0barSum_5Mev_Bin3", "D0bar invariant mass Bin3 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
+ fhInvMassD0barSum_5Mev_Bin3->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0barSum_5Mev_Bin3->SetMinimum(0);
+ fOutput3->Add(fhInvMassD0barSum_5Mev_Bin3);
+
+ fhInvMassSumAll_5Mev_Bin3 = new TH1F("fhInvMassSumAll_5Mev_Bin3", "D0/D0bar invariant mass Bin3 (good hyps); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
+ fhInvMassSumAll_5Mev_Bin3->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassSumAll_5Mev_Bin3->SetMinimum(0);
+ fOutput3->Add(fhInvMassSumAll_5Mev_Bin3);
+
+ fhInvMassD0Sum_10Mev_Bin4 = new TH1F("fhInvMassD0Sum_10Mev_Bin4", "D0 invariant mass Bin4 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
+ fhInvMassD0Sum_10Mev_Bin4->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0Sum_10Mev_Bin4->SetMinimum(0);
+ fOutput4->Add(fhInvMassD0Sum_10Mev_Bin4);
+
+ fhInvMassD0barSum_10Mev_Bin4 = new TH1F("fhInvMassD0barSum_10Mev_Bin4", "D0bar invariant mass Bin4 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
+ fhInvMassD0barSum_10Mev_Bin4->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0barSum_10Mev_Bin4->SetMinimum(0);
+ fOutput4->Add(fhInvMassD0barSum_10Mev_Bin4);
+
+ fhInvMassSumAll_10Mev_Bin4 = new TH1F("fhInvMassSumAll_10Mev_Bin4", "D0/D0bar invariant mass Bin4 (good hyps); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
+ fhInvMassSumAll_10Mev_Bin4->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassSumAll_10Mev_Bin4->SetMinimum(0);
+ fOutput4->Add(fhInvMassSumAll_10Mev_Bin4);
+
+ fhInvMassD0Sum_5Mev_Bin4 = new TH1F("fhInvMassD0Sum_5Mev_Bin4", "D0 invariant mass Bin4 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
+ fhInvMassD0Sum_5Mev_Bin4->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0Sum_5Mev_Bin4->SetMinimum(0);
+ fOutput4->Add(fhInvMassD0Sum_5Mev_Bin4);
+
+ fhInvMassD0barSum_5Mev_Bin4 = new TH1F("fhInvMassD0barSum_5Mev_Bin4", "D0bar invariant mass Bin4 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
+ fhInvMassD0barSum_5Mev_Bin4->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0barSum_5Mev_Bin4->SetMinimum(0);
+ fOutput4->Add(fhInvMassD0barSum_5Mev_Bin4);
+
+ fhInvMassSumAll_5Mev_Bin4 = new TH1F("fhInvMassSumAll_5Mev_Bin4", "D0/D0bar invariant mass Bin4 (good hyps); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
+ fhInvMassSumAll_5Mev_Bin4->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassSumAll_5Mev_Bin4->SetMinimum(0);
+ fOutput4->Add(fhInvMassSumAll_5Mev_Bin4);
+
+ fhInvMassD0Sum_10Mev_Bin5 = new TH1F("fhInvMassD0Sum_10Mev_Bin5", "D0 invariant mass Bin5 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
+ fhInvMassD0Sum_10Mev_Bin5->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0Sum_10Mev_Bin5->SetMinimum(0);
+ fOutput5->Add(fhInvMassD0Sum_10Mev_Bin5);
+
+ fhInvMassD0barSum_10Mev_Bin5 = new TH1F("fhInvMassD0barSum_10Mev_Bin5", "D0bar invariant mass Bin5 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
+ fhInvMassD0barSum_10Mev_Bin5->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0barSum_10Mev_Bin5->SetMinimum(0);
+ fOutput5->Add(fhInvMassD0barSum_10Mev_Bin5);
+
+ fhInvMassSumAll_10Mev_Bin5 = new TH1F("fhInvMassSumAll_10Mev_Bin5", "D0/D0bar invariant mass Bin5 (good hyps); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
+ fhInvMassSumAll_10Mev_Bin5->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassSumAll_10Mev_Bin5->SetMinimum(0);
+ fOutput5->Add(fhInvMassSumAll_10Mev_Bin5);
+
+ fhInvMassD0Sum_5Mev_Bin5 = new TH1F("fhInvMassD0Sum_5Mev_Bin5", "D0 invariant mass Bin5 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
+ fhInvMassD0Sum_5Mev_Bin5->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0Sum_5Mev_Bin5->SetMinimum(0);
+ fOutput5->Add(fhInvMassD0Sum_5Mev_Bin5);
+
+ fhInvMassD0barSum_5Mev_Bin5 = new TH1F("fhInvMassD0barSum_5Mev_Bin5", "D0bar invariant mass Bin5 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
+ fhInvMassD0barSum_5Mev_Bin5->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassD0barSum_5Mev_Bin5->SetMinimum(0);
+ fOutput5->Add(fhInvMassD0barSum_5Mev_Bin5);
+
+ fhInvMassSumAll_5Mev_Bin5 = new TH1F("fhInvMassSumAll_5Mev_Bin5", "D0/D0bar invariant mass Bin5 (good hyps); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
+ fhInvMassSumAll_5Mev_Bin5->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassSumAll_5Mev_Bin5->SetMinimum(0);
+ fOutput5->Add(fhInvMassSumAll_5Mev_Bin5);
+
+ fhInvMassMultipleOnly_Bin1 = new TH1F("fhInvMassMultipleOnly_Bin1", "D0/D0bar invariant mass Bin1 (good hyps) - Multple hyps accepted only; Inv. mass [GeV]; Entries/10 MeV",120,1.6,2.2);
+ fhInvMassMultipleOnly_Bin1->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassMultipleOnly_Bin1->SetMinimum(0);
+ fOutput->Add(fhInvMassMultipleOnly_Bin1);
+
+ fhInvMassMultipleOnly_Bin2 = new TH1F("fhInvMassMultipleOnly_Bin2", "D0/D0bar invariant mass Bin2 (good hyps) - Multple hyps accepted only; Inv. mass [GeV]; Entries/10 MeV",120,1.6,2.2);
+ fhInvMassMultipleOnly_Bin2->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassMultipleOnly_Bin2->SetMinimum(0);
+ fOutput2->Add(fhInvMassMultipleOnly_Bin2);
+
+ fhInvMassMultipleOnly_Bin3 = new TH1F("fhInvMassMultipleOnly_Bin3", "D0/D0bar invariant mass Bin3 (good hyps) - Multple hyps accepted only; Inv. mass [GeV]; Entries/10 MeV",120,1.6,2.2);
+ fhInvMassMultipleOnly_Bin3->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassMultipleOnly_Bin3->SetMinimum(0);
+ fOutput3->Add(fhInvMassMultipleOnly_Bin3);
+
+ fhInvMassMultipleOnly_Bin4 = new TH1F("fhInvMassMultipleOnly_Bin4", "D0/D0bar invariant mass Bin4 (good hyps) - Multple hyps accepted only; Inv. mass [GeV]; Entries/10 MeV",120,1.6,2.2);
+ fhInvMassMultipleOnly_Bin4->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassMultipleOnly_Bin4->SetMinimum(0);
+ fOutput4->Add(fhInvMassMultipleOnly_Bin4);
+
+ fhInvMassMultipleOnly_Bin5 = new TH1F("fhInvMassMultipleOnly_Bin5", "D0/D0bar invariant mass Bin5 (good hyps) - Multple hyps accepted only; Inv. mass [GeV]; Entries/10 MeV",120,1.6,2.2);
+ fhInvMassMultipleOnly_Bin5->Sumw2(); //Create structure to store sum of squares of weights
+ fhInvMassMultipleOnly_Bin5->SetMinimum(0);
+ fOutput5->Add(fhInvMassMultipleOnly_Bin5);
+
+ 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);
+ fOutput->Add(fScatterP4PID);
+
+ fPtVsY = new TH2F("fPtVsY", "Pt vs Y PPR Sel. Candidates; Pt [GeV/c]; Y",250,0.,25.,300,-3.,3.);
+ fPtVsY->SetMinimum(0);
+ fOutputC->Add(fPtVsY);
+
+ fPtVsYAll = new TH2F("fPtVsYAll", "Pt vs Y All Candidates; Pt [GeV/c]; Y",250,0.,25.,300,-3.,3.);
+ fPtVsYAll->SetMinimum(0);
+ fOutputC->Add(fPtVsYAll);
+
+ fEventCounter = new TH1F("fEventCounter", "N° of total events; NA; Events",1,0.,1.);
+ fEventCounter->SetMinimum(0);
+ fOutputC->Add(fEventCounter);
+
+ fCutDCA = new TH1F("fCutDCA", "DCA of candidate (couple); DCA [cm]; Entries/micron",500,0.,0.05);
+ fCutDCA->SetMinimum(0);
+ fOutputC->Add(fCutDCA);
+
+ fCutDCA3 = new TH1F("fCutDCA3", "DCA of candidate (trips); DCA [cm]; Entries/micron",500,0.,0.05);
+ fCutDCA3->SetMinimum(0);
+ fOutputC->Add(fCutDCA3);
+
+ fCutDCA2 = new TH1F("fCutDCA2", "DCA of candidate (quads1); DCA [cm]; Entries/micron",500,0.,0.05);
+ fCutDCA2->SetMinimum(0);
+ fOutputC->Add(fCutDCA2);
+
+ fCutDCA5 = new TH1F("fCutDCA5", "DCA of candidate (quads2); DCA [cm]; Entries/micron",500,0.,0.05);
+ fCutDCA5->SetMinimum(0);
+ fOutputC->Add(fCutDCA5);
+
+ fCutVertexDist2 = new TH1F("fCutVertexDist2", "Distance Vtx doubl.-Primary Vtx; Distance [cm]; Entries/15 micron",500,0.,0.75);
+ fCutVertexDist2->SetMinimum(0);
+ fOutputC->Add(fCutVertexDist2);
+
+ fCutVertexDist3 = new TH1F("fCutVertexDist3", "Distance Vtx trips-Primary Vtx; Distance [cm]; Entries/10 micron",500,0.,0.5);
+ fCutVertexDist3->SetMinimum(0);
+ fOutputC->Add(fCutVertexDist3);
+
+ fCutVertexDist4 = new TH1F("fCutVertexDist4", "Distance Vtx quads-Primary Vtx; Distance [cm]; Entries/5 micron",500,0.,0.25);
+ fCutVertexDist4->SetMinimum(0);
+ fOutputC->Add(fCutVertexDist4);
+
+ fCutCosinePoint = new TH1F("fCutCosinePoint", "Cosine of angle of pointing; Cos(Thetapt.); Entries/10^(-3)",250,0.75,1.);
+ fCutCosinePoint->SetMinimum(0);
+ fOutputC->Add(fCutCosinePoint);
+
+ fCutPt = new TH1F("fCutPt", "Pt of candidate D0; Pt [GeV/c]; Entries/5 MeV",3000,0.,15.);
+ fCutPt->SetMinimum(0);
+ fOutputC->Add(fCutPt);
+
+ fCutY = new TH1F("fCutY", "Y of candidate D0; Pt [GeV/c]; Entries/5 MeV",900,-9.,9.);
+ fCutY->SetMinimum(0);
+ fOutputC->Add(fCutY);
+
+ fPIDSel = new TH1F("fPIDSel", "Ratio of D0 selected by PID for Correction",3,0.,3.);
+ fPIDSel->SetMinimum(0);
+ fPIDSel->GetXaxis()->SetBinLabel(1,"D0allhyp All");
+ fPIDSel->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
+ fPIDSel->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
+
+ fPIDSel_Bin1 = new TH1F("fPIDSel_Bin1", "Ratio of D0 selected by PID for Correction",3,0.,3.);
+ fPIDSel_Bin1->SetMinimum(0);
+ fPIDSel_Bin1->GetXaxis()->SetBinLabel(1,"D0allhyp All");
+ fPIDSel_Bin1->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
+ fPIDSel_Bin1->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
+
+ fPIDSel_Bin2 = new TH1F("fPIDSel_Bin1", "Ratio of D0 selected by PID for Correction",3,0.,3.);
+ fPIDSel_Bin2->SetMinimum(0);
+ fPIDSel_Bin2->GetXaxis()->SetBinLabel(1,"D0allhyp All");
+ fPIDSel_Bin2->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
+ fPIDSel_Bin2->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
+
+ fPIDSel_Bin3 = new TH1F("fPIDSel_Bin1", "Ratio of D0 selected by PID for Correction",3,0.,3.);
+ fPIDSel_Bin3->SetMinimum(0);
+ fPIDSel_Bin3->GetXaxis()->SetBinLabel(1,"D0allhyp All");
+ fPIDSel_Bin3->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
+ fPIDSel_Bin3->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
+
+ fPIDSel_Bin4 = new TH1F("fPIDSel_Bin1", "Ratio of D0 selected by PID for Correction",3,0.,3.);
+ fPIDSel_Bin4->SetMinimum(0);
+ fPIDSel_Bin4->GetXaxis()->SetBinLabel(1,"D0allhyp All");
+ fPIDSel_Bin4->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
+ fPIDSel_Bin4->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
+
+ fPIDSel_Bin5 = new TH1F("fPIDSel_Bin1", "Ratio of D0 selected by PID for Correction",3,0.,3.);
+ fPIDSel_Bin5->SetMinimum(0);
+ fPIDSel_Bin5->GetXaxis()->SetBinLabel(1,"D0allhyp All");
+ fPIDSel_Bin5->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
+ fPIDSel_Bin5->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
+
+ fMultipleHyps = new TH1F("fMultipleHyps", "N. of hyp. accepted for each candidate (accounted N. times)",8,0.,8.);
+ fMultipleHyps->SetMinimum(0);
+ fMultipleHyps->GetXaxis()->SetBinLabel(1,"1 (PPR)");
+ fMultipleHyps->GetXaxis()->SetBinLabel(2,"2 (PPR)");
+ fMultipleHyps->GetXaxis()->SetBinLabel(3,"3 (PPR)");
+ fMultipleHyps->GetXaxis()->SetBinLabel(4,"4 (PPR)");
+ fMultipleHyps->GetXaxis()->SetBinLabel(5,"1 (PPR+PID)");
+ fMultipleHyps->GetXaxis()->SetBinLabel(6,"2 (PPR+PID)");
+ fMultipleHyps->GetXaxis()->SetBinLabel(7,"3 (PPR+PID)");
+ fMultipleHyps->GetXaxis()->SetBinLabel(8,"4 (PPR+PID)");
+
+ fMultipleHypsType = new TH1F("fMultipleHypsType", "Type of hyp. accepted for each candidate",8,0.,8.);
+ fMultipleHypsType->SetMinimum(0);
+ fMultipleHypsType->GetXaxis()->SetBinLabel(1,"D0");
+ fMultipleHypsType->GetXaxis()->SetBinLabel(2,"D0bar");
+ fMultipleHypsType->GetXaxis()->SetBinLabel(3,"2D0");
+ fMultipleHypsType->GetXaxis()->SetBinLabel(4,"2D0bar");
+ fMultipleHypsType->GetXaxis()->SetBinLabel(5,"D0+D0bar");
+ fMultipleHypsType->GetXaxis()->SetBinLabel(6,"2D0+D0bar");
+ fMultipleHypsType->GetXaxis()->SetBinLabel(7,"D0+2D0bar");
+ fMultipleHypsType->GetXaxis()->SetBinLabel(8,"2D0+2D0bar");
+
+ fOutputC->Add(fMultipleHyps);
+ fOutputC->Add(fMultipleHypsType);
+
+ fOutputC->Add(fPIDSel);
+ fOutput->Add(fPIDSel_Bin1);
+ fOutput2->Add(fPIDSel_Bin2);
+ fOutput3->Add(fPIDSel_Bin3);
+ fOutput4->Add(fPIDSel_Bin4);
+ fOutput5->Add(fPIDSel_Bin5);
+
+ fPtSel = new TH1F("fPtSel", "Pt of candidates accepted; Pt [GeV/c]; Entries/10 MeV",2000,0.,20.);
+ fPtSel->SetMinimum(0);
+ fOutputC->Add(fPtSel);
+
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSESelectHF4Prong::UserExec(Option_t */*option*/)
+{
+ // Execute analysis for current event:
+ // heavy flavor candidates selection and histograms
+
+ AliAODEvent *aodIn = dynamic_cast<AliAODEvent*> (InputEvent());
+
+ TClonesArray *inputArrayCharm4Prong = 0;
+
+ if(!aodIn && AODEvent() && IsStandardAOD()) {
+ // In case there is an AOD handler writing a standard AOD, use the AOD
+ // event in memory rather than the input (ESD) event.
+ aodIn = dynamic_cast<AliAODEvent*> (AODEvent());
+ // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
+ // have to taken from the AOD event hold by the AliAODExtension
+ AliAODHandler* aodHandler = (AliAODHandler*)
+ ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
+ if(aodHandler->GetExtensions()) {
+ AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
+ AliAODEvent *aodFromExt = ext->GetAOD();
+ // load D0 candidates
+ inputArrayCharm4Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
+ }
+ } else {
+ // load D0 candidates
+ inputArrayCharm4Prong=(TClonesArray*)aodIn->GetList()->FindObject("Charm4Prong");
+ }
+
+ if(!inputArrayCharm4Prong) {
+ printf("AliAnalysisTaskSESelectHF4Prong::UserExec: D0to3Kpi branch not found!\n");
+ return;
+ }
+
+ //print event info
+// aodIn->GetHeader()->Print();
+
+ //Event counter ++
+ fEventCounter->Fill(0);
+
+ // fix for temporary bug in ESDfilter
+ // the AODs with null vertex pointer didn't pass the PhysSel
+ if(!aodIn->GetPrimaryVertex()) return;
+
+ // primary vertex
+ AliAODVertex *vtx1 = (AliAODVertex*)aodIn->GetPrimaryVertex();
+// vtx1->Print();
+
+ // make trkIDtoEntry register (temporary)
+ Int_t trkIDtoEntry[100000];
+ for(Int_t it=0;it<aodIn->GetNumberOfTracks();it++) {
+ AliAODTrack *track = aodIn->GetTrack(it);
+ trkIDtoEntry[track->GetID()]=it;
+ }
+
+ Int_t iOutVerticesHF=0,iOutCharm4Prong=0;
+ fVerticesHFTClArr->Delete();
+ iOutVerticesHF = fVerticesHFTClArr->GetEntriesFast();
+ TClonesArray &verticesHFRef = *fVerticesHFTClArr;
+ fCharm4ProngTClArr->Delete();
+ iOutCharm4Prong = fCharm4ProngTClArr->GetEntriesFast();
+ TClonesArray &aodCharm4ProngRef = *fCharm4ProngTClArr;
+
+ // loop over D0->K3pi candidates
+ Int_t nInCharm4Prong = inputArrayCharm4Prong->GetEntriesFast();
+ printf("Number of D0->K3pi: %d\n",nInCharm4Prong);
+
+ for (Int_t iCharm4Prong = 0; iCharm4Prong < nInCharm4Prong; iCharm4Prong++) {
+ AliAODRecoDecayHF4Prong *dIn = (AliAODRecoDecayHF4Prong*)inputArrayCharm4Prong->UncheckedAt(iCharm4Prong);
+ Bool_t unsetvtx=kFALSE;
+
+ if(!dIn->GetOwnPrimaryVtx()) {
+ dIn->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
+ unsetvtx=kTRUE;
+ }
+
+ //fill histos of cuts
+ Double_t dca = dIn->GetDCA();
+ Double_t dist2 = dIn->GetDist12toPrim();
+ Double_t dist3 = dIn->GetDist3toPrim();
+ Double_t dist4 = dIn->GetDist4toPrim();
+ Double_t cosine = dIn->CosPointingAngle();
+ Double_t ptPart = dIn->Pt();
+ Double_t yPart = dIn->YD0();
+ Double_t dcatrip = dIn->GetDCA(3);
+ 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);
+ fCutDCA5->Fill(dcaquad2);
+ fCutVertexDist2->Fill(dist2);
+ fCutVertexDist3->Fill(dist3);
+ fCutVertexDist4->Fill(dist4);
+ fCutCosinePoint->Fill(cosine);
+ fCutPt->Fill(ptPart);
+ fCutY->Fill(yPart);
+ fPtVsYAll->Fill(ptPart,yPart);
+
+ //flags initialization
+ fSelected = fCuts->IsSelected(dIn,AliRDHFCuts::kCandidate);
+ Int_t selD01 = fCuts->D01Selected(dIn,AliRDHFCuts::kCandidate);
+ Int_t selD02 = fCuts->D02Selected(dIn,AliRDHFCuts::kCandidate);
+ Int_t selD0bar1 = fCuts->D0bar1Selected(dIn,AliRDHFCuts::kCandidate);
+ Int_t selD0bar2 = fCuts->D0bar2Selected(dIn,AliRDHFCuts::kCandidate);
+ Int_t flagAccLim = 1;
+
+ //Limited Acceptance
+ if(ptPart > 5.) {
+ if (TMath::Abs(yPart) > 0.8) flagAccLim = 0;
+ }
+ else {
+ Double_t maxFiducialY = -0.2/15*ptPart*ptPart+1.9/15*ptPart+0.5;
+ Double_t minFiducialY = 0.2/15*ptPart*ptPart-1.9/15*ptPart-0.5;
+ if (yPart < minFiducialY || yPart > maxFiducialY) flagAccLim = 0;;
+ }
+
+ //number of CANDIDATES (regardless of hypotheses) passing PPR
+ if(fSelected==1||fSelected==2||fSelected==3) {
+ fPtVsY->Fill(ptPart,yPart);
+ fPIDSel->Fill(0);
+ if (ptPart >= ptBinH[0] && ptPart < ptBinH[1]) fPIDSel_Bin1->Fill(0);
+ else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fPIDSel_Bin2->Fill(0);
+ else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fPIDSel_Bin3->Fill(0);
+ else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fPIDSel_Bin4->Fill(0);
+ else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fPIDSel_Bin5->Fill(0);
+ }
+
+
+ PostData(6,fOutputC);
+
+ //selection
+ if((fSelected==1||fSelected==2||fSelected==3) && flagAccLim == 1) {
+ // get daughter AOD tracks
+ AliAODTrack *trk0 = (AliAODTrack*)dIn->GetDaughter(0);
+ AliAODTrack *trk1 = (AliAODTrack*)dIn->GetDaughter(1);
+ AliAODTrack *trk2 = (AliAODTrack*)dIn->GetDaughter(2);
+ AliAODTrack *trk3 = (AliAODTrack*)dIn->GetDaughter(3);
+ if(!trk0 || !trk1 || !trk2 || !trk3) {
+ trk0=aodIn->GetTrack(trkIDtoEntry[dIn->GetProngID(0)]);
+ trk1=aodIn->GetTrack(trkIDtoEntry[dIn->GetProngID(1)]);
+ trk2=aodIn->GetTrack(trkIDtoEntry[dIn->GetProngID(2)]);
+ trk3=aodIn->GetTrack(trkIDtoEntry[dIn->GetProngID(3)]);
+ }
+ printf("pt of positive track #1: %f\n",trk0->Pt());
+ printf("pt of negative track #1: %f\n",trk1->Pt());
+ printf("pt of positive track #2: %f\n",trk2->Pt());
+ printf("pt of negative track #2: %f\n",trk3->Pt());
+
+ dIn->InvMassD0(fmassD0);
+ dIn->InvMassD0bar(fmassD0bar);
+
+ //fill histos (combining selection form cuts & PID (with rho information))
+ 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);
+
+ //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]) fPIDSel_Bin1->Fill(1);
+ else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fPIDSel_Bin2->Fill(1);
+ else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fPIDSel_Bin3->Fill(1);
+ else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fPIDSel_Bin4->Fill(1);
+ else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fPIDSel_Bin5->Fill(1);
+ }
+
+ //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);
+
+ //combine PPD + PID cuts
+ if (selD01 == 1 && pid1==1) hypD01 = 1;
+ if (selD02 == 1 && pid2==1) hypD02 = 1;
+ if (selD0bar1 == 1 && pidbar1==1) hypD0bar1 = 1;
+ if (selD0bar2 == 1 && pidbar2==1) hypD0bar2 = 1;
+
+ //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]) fPIDSel_Bin1->Fill(2);
+ else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fPIDSel_Bin2->Fill(2);
+ else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fPIDSel_Bin3->Fill(2);
+ else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fPIDSel_Bin4->Fill(2);
+ else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fPIDSel_Bin5->Fill(2);
+ }
+
+ //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);
+
+ //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]) fhInvMassMultipleOnly_Bin1->Fill(fmassD0[0]);
+ else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fhInvMassMultipleOnly_Bin2->Fill(fmassD0[0]);
+ else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fhInvMassMultipleOnly_Bin3->Fill(fmassD0[0]);
+ else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fhInvMassMultipleOnly_Bin4->Fill(fmassD0[0]);
+ else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fhInvMassMultipleOnly_Bin5->Fill(fmassD0[0]);
+ }
+ // D02 hyp.
+ if(hypD02==1) {
+ if (ptPart < ptBinH[1]) fhInvMassMultipleOnly_Bin1->Fill(fmassD0[1]);
+ else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fhInvMassMultipleOnly_Bin2->Fill(fmassD0[1]);
+ else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fhInvMassMultipleOnly_Bin3->Fill(fmassD0[1]);
+ else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fhInvMassMultipleOnly_Bin4->Fill(fmassD0[1]);
+ else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fhInvMassMultipleOnly_Bin5->Fill(fmassD0[1]);
+ }
+ // D0bar1 hyp.
+ if(hypD0bar1==1) {
+ if (ptPart < ptBinH[1]) fhInvMassMultipleOnly_Bin1->Fill(fmassD0bar[0]);
+ else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fhInvMassMultipleOnly_Bin2->Fill(fmassD0bar[0]);
+ else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fhInvMassMultipleOnly_Bin3->Fill(fmassD0bar[0]);
+ else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fhInvMassMultipleOnly_Bin4->Fill(fmassD0bar[0]);
+ else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fhInvMassMultipleOnly_Bin5->Fill(fmassD0bar[0]);
+ }
+ // D0bar2 hyp.
+ if(hypD0bar2==1) {
+ if (ptPart < ptBinH[1]) fhInvMassMultipleOnly_Bin1->Fill(fmassD0bar[1]);
+ else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fhInvMassMultipleOnly_Bin2->Fill(fmassD0bar[1]);
+ else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fhInvMassMultipleOnly_Bin3->Fill(fmassD0bar[1]);
+ else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fhInvMassMultipleOnly_Bin4->Fill(fmassD0bar[1]);
+ else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fhInvMassMultipleOnly_Bin5->Fill(fmassD0bar[1]);
+ }
+ }
+
+ }
+//
+//end of auxiliary section
+//
+
+
+ //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]) {
+
+ // D01 hyp.
+ if(hypD01==1) {
+ fPtSel->Fill(ptPart);
+ fScatterP4PID->Fill(trk1->Pt(),trk3->Pt());
+ if (ptPart < ptBinH[1]) {fhInvMassD0Sum_10Mev_Bin1->Fill(fmassD0[0]);
+ fhInvMassD0Sum_5Mev_Bin1->Fill(fmassD0[0]);
+ fhInvMassSumAll_10Mev_Bin1->Fill(fmassD0[0]);
+ fhInvMassSumAll_5Mev_Bin1->Fill(fmassD0[0]);}
+ else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) {fhInvMassD0Sum_10Mev_Bin2->Fill(fmassD0[0]);
+ fhInvMassD0Sum_5Mev_Bin2->Fill(fmassD0[0]);
+ fhInvMassSumAll_10Mev_Bin2->Fill(fmassD0[0]);
+ fhInvMassSumAll_5Mev_Bin2->Fill(fmassD0[0]);}
+ else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) {fhInvMassD0Sum_10Mev_Bin3->Fill(fmassD0[0]);
+ fhInvMassD0Sum_5Mev_Bin3->Fill(fmassD0[0]);
+ fhInvMassSumAll_10Mev_Bin3->Fill(fmassD0[0]);
+ fhInvMassSumAll_5Mev_Bin3->Fill(fmassD0[0]);}
+ else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) {fhInvMassD0Sum_10Mev_Bin4->Fill(fmassD0[0]);
+ fhInvMassD0Sum_5Mev_Bin4->Fill(fmassD0[0]);
+ fhInvMassSumAll_10Mev_Bin4->Fill(fmassD0[0]);
+ fhInvMassSumAll_5Mev_Bin4->Fill(fmassD0[0]);}
+ else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) {fhInvMassD0Sum_10Mev_Bin5->Fill(fmassD0[0]);
+ fhInvMassD0Sum_5Mev_Bin5->Fill(fmassD0[0]);
+ fhInvMassSumAll_10Mev_Bin5->Fill(fmassD0[0]);
+ fhInvMassSumAll_5Mev_Bin5->Fill(fmassD0[0]);}
+ }
+ // D02 hyp.
+ if(hypD02==1) {
+ fPtSel->Fill(ptPart);
+ fScatterP4PID->Fill(trk3->Pt(),trk1->Pt());
+ if (ptPart < ptBinH[1]) {fhInvMassD0Sum_10Mev_Bin1->Fill(fmassD0[1]);
+ fhInvMassD0Sum_5Mev_Bin1->Fill(fmassD0[1]);
+ fhInvMassSumAll_10Mev_Bin1->Fill(fmassD0[1]);
+ fhInvMassSumAll_5Mev_Bin1->Fill(fmassD0[1]);}
+ else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) {fhInvMassD0Sum_10Mev_Bin2->Fill(fmassD0[1]);
+ fhInvMassD0Sum_5Mev_Bin2->Fill(fmassD0[1]);
+ fhInvMassSumAll_10Mev_Bin2->Fill(fmassD0[1]);
+ fhInvMassSumAll_5Mev_Bin2->Fill(fmassD0[1]);}
+ else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) {fhInvMassD0Sum_10Mev_Bin3->Fill(fmassD0[1]);
+ fhInvMassD0Sum_5Mev_Bin3->Fill(fmassD0[1]);
+ fhInvMassSumAll_10Mev_Bin3->Fill(fmassD0[1]);
+ fhInvMassSumAll_5Mev_Bin3->Fill(fmassD0[1]);}
+ else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) {fhInvMassD0Sum_10Mev_Bin4->Fill(fmassD0[1]);
+ fhInvMassD0Sum_5Mev_Bin4->Fill(fmassD0[1]);
+ fhInvMassSumAll_10Mev_Bin4->Fill(fmassD0[1]);
+ fhInvMassSumAll_5Mev_Bin4->Fill(fmassD0[1]);}
+ else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) {fhInvMassD0Sum_10Mev_Bin5->Fill(fmassD0[1]);
+ fhInvMassD0Sum_5Mev_Bin5->Fill(fmassD0[1]);
+ fhInvMassSumAll_10Mev_Bin5->Fill(fmassD0[1]);
+ fhInvMassSumAll_5Mev_Bin5->Fill(fmassD0[1]);}
+ }
+ // D0bar1 hyp.
+ if(hypD0bar1==1) {
+ fPtSel->Fill(ptPart);
+ fScatterP4PID->Fill(trk0->Pt(),trk2->Pt());
+ if (ptPart < ptBinH[1]) {fhInvMassD0barSum_10Mev_Bin1->Fill(fmassD0bar[0]);
+ fhInvMassD0barSum_5Mev_Bin1->Fill(fmassD0bar[0]);
+ fhInvMassSumAll_10Mev_Bin1->Fill(fmassD0bar[0]);
+ fhInvMassSumAll_5Mev_Bin1->Fill(fmassD0bar[0]);}
+ else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) {fhInvMassD0barSum_10Mev_Bin2->Fill(fmassD0bar[0]);
+ fhInvMassD0barSum_5Mev_Bin2->Fill(fmassD0bar[0]);
+ fhInvMassSumAll_10Mev_Bin2->Fill(fmassD0bar[0]);
+ fhInvMassSumAll_5Mev_Bin2->Fill(fmassD0bar[0]);}
+ else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) {fhInvMassD0barSum_10Mev_Bin3->Fill(fmassD0bar[0]);
+ fhInvMassD0barSum_5Mev_Bin3->Fill(fmassD0bar[0]);
+ fhInvMassSumAll_10Mev_Bin3->Fill(fmassD0bar[0]);
+ fhInvMassSumAll_5Mev_Bin3->Fill(fmassD0bar[0]);}
+ else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) {fhInvMassD0barSum_10Mev_Bin4->Fill(fmassD0bar[0]);
+ fhInvMassD0barSum_5Mev_Bin4->Fill(fmassD0bar[0]);
+ fhInvMassSumAll_10Mev_Bin4->Fill(fmassD0bar[0]);
+ fhInvMassSumAll_5Mev_Bin4->Fill(fmassD0bar[0]);}
+ else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) {fhInvMassD0barSum_10Mev_Bin5->Fill(fmassD0bar[0]);
+ fhInvMassD0barSum_5Mev_Bin5->Fill(fmassD0bar[0]);
+ fhInvMassSumAll_10Mev_Bin5->Fill(fmassD0bar[0]);
+ fhInvMassSumAll_5Mev_Bin5->Fill(fmassD0bar[0]);}
+ }
+ // D0bar2 hyp.
+ if(hypD0bar2==1) {
+ fPtSel->Fill(ptPart);
+ fScatterP4PID->Fill(trk2->Pt(),trk0->Pt());
+ if (ptPart < ptBinH[1]) {fhInvMassD0barSum_10Mev_Bin1->Fill(fmassD0bar[1]);
+ fhInvMassD0barSum_5Mev_Bin1->Fill(fmassD0bar[1]);
+ fhInvMassSumAll_10Mev_Bin1->Fill(fmassD0bar[1]);
+ fhInvMassSumAll_5Mev_Bin1->Fill(fmassD0bar[1]);}
+ else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) {fhInvMassD0barSum_10Mev_Bin2->Fill(fmassD0bar[1]);
+ fhInvMassD0barSum_5Mev_Bin2->Fill(fmassD0bar[1]);
+ fhInvMassSumAll_10Mev_Bin2->Fill(fmassD0bar[1]);
+ fhInvMassSumAll_5Mev_Bin2->Fill(fmassD0bar[1]);}
+ else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) {fhInvMassD0barSum_10Mev_Bin3->Fill(fmassD0bar[1]);
+ fhInvMassD0barSum_5Mev_Bin3->Fill(fmassD0bar[1]);
+ fhInvMassSumAll_10Mev_Bin3->Fill(fmassD0bar[1]);
+ fhInvMassSumAll_5Mev_Bin3->Fill(fmassD0bar[1]);}
+ else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) {fhInvMassD0barSum_10Mev_Bin4->Fill(fmassD0bar[1]);
+ fhInvMassD0barSum_5Mev_Bin4->Fill(fmassD0bar[1]);
+ fhInvMassSumAll_10Mev_Bin4->Fill(fmassD0bar[1]);
+ fhInvMassSumAll_5Mev_Bin4->Fill(fmassD0bar[1]);}
+ else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) {fhInvMassD0barSum_10Mev_Bin5->Fill(fmassD0bar[1]);
+ fhInvMassD0barSum_5Mev_Bin5->Fill(fmassD0bar[1]);
+ fhInvMassSumAll_10Mev_Bin5->Fill(fmassD0bar[1]);
+ fhInvMassSumAll_5Mev_Bin5->Fill(fmassD0bar[1]);}
+ }
+ }
+
+ PostData(1,fOutput);
+ PostData(2,fOutput2);
+ PostData(3,fOutput3);
+ PostData(4,fOutput4);
+ PostData(5,fOutput5);
+
+ // HERE ONE COULD RECALCULATE THE VERTEX USING THE KF PACKAGE
+
+ // clone candidate for output AOD
+ if(hypD01||hypD02||hypD0bar1||hypD0bar2) {
+ AliAODVertex *v = new(verticesHFRef[iOutVerticesHF++])
+ AliAODVertex(*(dIn->GetSecondaryVtx()));
+ AliAODRecoDecayHF4Prong *dOut=new(aodCharm4ProngRef[iOutCharm4Prong++])
+ AliAODRecoDecayHF4Prong(*dIn);
+ dOut->SetSecondaryVtx(v);
+ 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
+
+ printf("Number of selected D0->K3pi: %d\n",iOutCharm4Prong);
+
+ return;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSESelectHF4Prong::Terminate(Option_t */*option*/)
+{
+ // Terminate analysis
+ //
+ if(fDebug > 1) printf("AnalysisTaskSESelectHF4Prong: Terminate() \n");
+
+ fOutput = dynamic_cast<TList*> (GetOutputData(1));
+ if (!fOutput) {
+ printf("ERROR: fOutput not available\n");
+ return;
+ }
+ fOutput2 = dynamic_cast<TList*> (GetOutputData(2));
+ if (!fOutput2) {
+ printf("ERROR: fOutput not available\n");
+ return;
+ }
+ fOutput3 = dynamic_cast<TList*> (GetOutputData(3));
+ if (!fOutput3) {
+ printf("ERROR: fOutput not available\n");
+ return;
+ }
+ fOutput4 = dynamic_cast<TList*> (GetOutputData(4));
+ if (!fOutput4) {
+ printf("ERROR: fOutput not available\n");
+ return;
+ }
+ if (!fOutput5) {
+ printf("ERROR: fOutput not available\n");
+ return;
+ }
+ fOutputC = dynamic_cast<TList*> (GetOutputData(5));
+ if (!fOutputC) {
+ printf("ERROR: fOutputC not available\n");
+ return;
+ }
+
+ fhInvMassD0Sum_10Mev_Bin1 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_10Mev_Bin1"));
+ fhInvMassD0barSum_10Mev_Bin1 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_10Mev_Bin1"));
+ fhInvMassSumAll_10Mev_Bin1 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_10Mev_Bin1"));
+ fhInvMassD0Sum_5Mev_Bin1 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_5Mev_Bin1"));
+ fhInvMassD0barSum_5Mev_Bin1 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_5Mev_Bin1"));
+ fhInvMassSumAll_5Mev_Bin1 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_5Mev_Bin1"));
+ fhInvMassD0Sum_10Mev_Bin2 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_10Mev_Bin2"));
+ fhInvMassD0barSum_10Mev_Bin2 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_10Mev_Bin2"));
+ fhInvMassSumAll_10Mev_Bin2 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_10Mev_Bin2"));
+ fhInvMassD0Sum_5Mev_Bin2 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_5Mev_Bin2"));
+ fhInvMassD0barSum_5Mev_Bin2 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_5Mev_Bin2"));
+ fhInvMassSumAll_5Mev_Bin2 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_5Mev_Bin2"));
+ fhInvMassD0Sum_10Mev_Bin3 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_10Mev_Bin3"));
+ fhInvMassD0barSum_10Mev_Bin3 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_10Mev_Bin3"));
+ fhInvMassSumAll_10Mev_Bin3 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_10Mev_Bin3"));
+ fhInvMassD0Sum_5Mev_Bin3 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_5Mev_Bin3"));
+ fhInvMassD0barSum_5Mev_Bin3 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_5Mev_Bin3"));
+ fhInvMassSumAll_5Mev_Bin3 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_5Mev_Bin3"));
+ fhInvMassD0Sum_10Mev_Bin4 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_10Mev_Bin4"));
+ fhInvMassD0barSum_10Mev_Bin4 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_10Mev_Bin4"));
+ fhInvMassSumAll_10Mev_Bin4 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_10Mev_Bin4"));
+ fhInvMassD0Sum_5Mev_Bin4 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_5Mev_Bin4"));
+ fhInvMassD0barSum_5Mev_Bin4 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_5Mev_Bin4"));
+ fhInvMassSumAll_5Mev_Bin4 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_5Mev_Bin4"));
+ fhInvMassD0Sum_10Mev_Bin5 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_10Mev_Bin5"));
+ fhInvMassD0barSum_10Mev_Bin5 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_10Mev_Bin5"));
+ fhInvMassSumAll_10Mev_Bin5 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_10Mev_Bin5"));
+ fhInvMassD0Sum_5Mev_Bin5 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_5Mev_Bin5"));
+ fhInvMassD0barSum_5Mev_Bin5 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_5Mev_Bin5"));
+ fhInvMassSumAll_5Mev_Bin5 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_5Mev_Bin5"));
+
+ fScatterP4PID = dynamic_cast<TH2F*>(fOutput->FindObject("fScatterP4PID"));
+ fPtVsY = dynamic_cast<TH2F*>(fOutputC->FindObject("fPtVsY"));
+ fPtVsYAll = dynamic_cast<TH2F*>(fOutputC->FindObject("fPtVsYAll"));
+
+ fEventCounter = dynamic_cast<TH1F*>(fOutputC->FindObject("fEventCounter"));
+
+ fCutDCA = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutDCA"));
+ fCutDCA3 = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutDCA3"));
+ fCutDCA2 = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutDCA2"));
+ fCutDCA5 = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutDCA5"));
+ fCutVertexDist2 = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutVertexDist2"));
+ fCutVertexDist3 = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutVertexDist3"));
+ fCutVertexDist4 = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutVertexDist4"));
+ fCutCosinePoint = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutCosinePoint"));
+
+ fCutPt = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutPt"));
+ fCutY = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutY"));
+ fPIDSel = dynamic_cast<TH1F*>(fOutputC->FindObject("fPIDSel"));
+ fPIDSel_Bin1 = dynamic_cast<TH1F*>(fOutputC->FindObject("fPIDSel_Bin1"));
+ fPIDSel_Bin2 = dynamic_cast<TH1F*>(fOutputC->FindObject("fPIDSel_Bin2"));
+ fPIDSel_Bin3 = dynamic_cast<TH1F*>(fOutputC->FindObject("fPIDSel_Bin3"));
+ fPIDSel_Bin4 = dynamic_cast<TH1F*>(fOutputC->FindObject("fPIDSel_Bin4"));
+ fPIDSel_Bin5 = dynamic_cast<TH1F*>(fOutputC->FindObject("fPIDSel_Bin5"));
+}
+
* provided "as is" without express or implied warranty. *\r
**************************************************************************/\r
\r
-/* $Id$ */\r
-\r
/////////////////////////////////////////////////////////////\r
//\r
// Class for cuts on AOD reconstructed D0->Kpipipi\r
//\r
-// Author: r.romita@gsi.de, andrea.dainese@pd.infn.it\r
+// Author: r.romita@gsi.de, andrea.dainese@pd.infn.it,
+// fabio.colamaria@ba.infn.it\r
/////////////////////////////////////////////////////////////\r
\r
#include <TDatabasePDG.h>\r
#include "AliAODRecoDecayHF4Prong.h"\r
#include "AliAODTrack.h"\r
#include "AliESDtrack.h"\r
+#include "AliAODPidHF.h"
\r
ClassImp(AliRDHFCutsD0toKpipipi)\r
\r
\r
Double_t ptD=d->Pt();\r
if(ptD<fMinPtCand) return 0;\r
- if(ptD>fMaxPtCand) return 0;\r
+ if(ptD>fMaxPtCand) return 0;
\r
// selection on daughter tracks \r
if(selectionLevel==AliRDHFCuts::kAll || \r
selectionLevel==AliRDHFCuts::kCandidate) {\r
\r
Int_t ptbin=PtBin(d->Pt());\r
- \r
- Int_t okD0=1,okD0bar=1; \r
+
+ Int_t okD0=1,okD0bar=1; \r
Double_t mD0[2],mD0bar[2];\r
Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
\r
TMath::Abs(mD0bar[1]-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0bar = 0;\r
if(!okD0 && !okD0bar) return 0;\r
\r
- if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) return 0;\r
+ if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]
+ || d->GetDCA(3) > fCutsRD[GetGlobalIndex(1,ptbin)]
+ || d->GetDCA(2) > fCutsRD[GetGlobalIndex(1,ptbin)]
+ || d->GetDCA(5) > fCutsRD[GetGlobalIndex(1,ptbin)]) return 0;\r
if(d->GetDist12toPrim() < fCutsRD[GetGlobalIndex(2,ptbin)]) return 0;\r
if(d->GetDist3toPrim() < fCutsRD[GetGlobalIndex(3,ptbin)]) return 0;\r
if(d->GetDist4toPrim() < fCutsRD[GetGlobalIndex(4,ptbin)]) return 0;\r
if(d->CosPointingAngle() < fCutsRD[GetGlobalIndex(5,ptbin)]) return 0;\r
if(d->Pt() < fCutsRD[GetGlobalIndex(6,ptbin)]) return 0;\r
if(!d->CutRhoMass(mD0,mD0bar,fCutsRD[GetGlobalIndex(0,ptbin)],fCutsRD[GetGlobalIndex(7,ptbin)])) return 0;\r
-\r
+
if (okD0) returnvalue=1; //cuts passed as D0\r
if (okD0bar) returnvalue=2; //cuts passed as D0bar\r
if (okD0 && okD0bar) returnvalue=3; //cuts passed as D0 and D0bar\r
}\r
-\r
+
+ // selection on PID (from AliAODPidHF)\r
+ if(selectionLevel==AliRDHFCuts::kAll || \r
+ selectionLevel==AliRDHFCuts::kPID) {\r
+
+ Int_t selD01 = D01Selected(d,AliRDHFCuts::kCandidate);
+ Int_t selD02 = D02Selected(d,AliRDHFCuts::kCandidate);
+ Int_t selD0bar1 = D0bar1Selected(d,AliRDHFCuts::kCandidate);
+ Int_t selD0bar2 = D0bar2Selected(d,AliRDHFCuts::kCandidate);
+
+ Int_t d01PID = 0, d02PID = 0, d0bar1PID = 0, d0bar2PID = 0;
+ returnvalue = IsSelectedFromPID(d, &d01PID, &d02PID, &d0bar1PID, &d0bar2PID); //This returnvalue is dummy! Now it's modified as it must be!
+
+returnvalue = 0;
+
+ if((selD01 == 1 && d01PID == 1)||(selD02 == 1 && d02PID == 1)||(selD0bar1 == 1 && d0bar1PID == 1)||(selD0bar2 == 1 && d0bar2PID == 1)) returnvalue = 1;\r
+ }
+
return returnvalue;\r
+}
+\r
+//---------------------------------------------------------------------------\r
+Int_t AliRDHFCutsD0toKpipipi::IsSelectedFromPID(AliAODRecoDecayHF4Prong *d, Int_t *hyp1, Int_t *hyp2, Int_t *hyp3, Int_t *hyp4) {\r
+ //\r
+ // Apply selection (using AliAODPidHF methods)
+ // Mass hypothesis true if each particle is at least compatible with specie of hypothesis\r
+ // \r
+
+ Int_t output=0;
+
+ Int_t matchK[4], matchPi[4];
+ Double_t ptlimit[2] = {0.6,0.8};\r
+ AliAODTrack* trk[4];
+ trk[0] = (AliAODTrack*)d->GetDaughter(0);
+ trk[1] = (AliAODTrack*)d->GetDaughter(1);
+ trk[2] = (AliAODTrack*)d->GetDaughter(2);
+ trk[3] = (AliAODTrack*)d->GetDaughter(3);
+// if(!trk[0] || !trk[1] || !trk[2] || !trk[3]) { //REMOVED (needs also AliAODEvent to be passed, here and in IsSelected method)
+// trk[0]=aodIn->GetTrack(trkIDtoEntry[d->GetProngID(0)]);
+// trk[1]=aodIn->GetTrack(trkIDtoEntry[d->GetProngID(1)]);
+// trk[2]=aodIn->GetTrack(trkIDtoEntry[d->GetProngID(2)]);
+// trk[3]=aodIn->GetTrack(trkIDtoEntry[d->GetProngID(3)]);}
+
+ AliAODPidHF* PidObj = new AliAODPidHF();
+ PidObj->SetAsym(kTRUE);
+ PidObj->SetPLimit(ptlimit);
+ PidObj->SetSigma(0,2.); //TPC sigma, in three pT ranges
+ PidObj->SetSigma(1,1.);
+ PidObj->SetSigma(2,0.);
+ PidObj->SetSigma(3,2.); //TOF sigma, whole pT range
+ PidObj->SetTPC(kTRUE);
+ PidObj->SetTOF(kTRUE);
+
+ for(Int_t ii=0; ii<4; ii++) {
+ PidObj->SetSigma(0,2.);
+ matchK[ii] = PidObj->MatchTPCTOF(trk[ii],1,3,kTRUE); //arguments: track, mode, particle#, compatibility allowed
+ PidObj->SetSigma(0,2.);
+ matchPi[ii] = PidObj->MatchTPCTOF(trk[ii],1,2,kTRUE);
+ }
+
+ //Rho invariant mass under various hypotheses (to be matched with PID infos in order to accet the candidate)
+ Int_t d01rho03 = 0, d01rho23 = 0, d02rho01 = 0, d02rho12 = 0, d0bar1rho12 = 0, d0bar1rho23 = 0, d0bar2rho01 = 0, d0bar2rho03 = 0;
+ if(TMath::Abs(0.775 - d->InvMassRho(0,3))<0.1) {d01rho03 = 1; d0bar2rho03 = 1;}
+ if(TMath::Abs(0.775 - d->InvMassRho(2,3))<0.1) {d01rho23 = 1; d0bar1rho23 = 1;}
+ if(TMath::Abs(0.775 - d->InvMassRho(0,1))<0.1) {d02rho01 = 1; d0bar2rho01 = 1;}
+ if(TMath::Abs(0.775 - d->InvMassRho(1,2))<0.1) {d02rho12 = 1; d0bar1rho12 = 1;}
+ Int_t d01rho = 0, d02rho = 0, d0bar1rho = 0, d0bar2rho = 0;
+ if(d01rho03==1||d01rho23==1) d01rho = 1;
+ if(d02rho01==1||d02rho12==1) d02rho = 1;
+ if(d0bar1rho12==1||d0bar1rho23==1) d0bar1rho = 1;
+ if(d0bar2rho01==1||d0bar2rho03==1) d0bar2rho = 1;
+
+ //This way because there could be multiple hypotheses accepted\r
+ if(d01rho==1 && (matchK[1]>=0 && matchPi[0]>=0 && matchPi[2]>=0 && matchPi[3]>=0)) {*hyp1 = 1; output = 1;} //d01 hyp
+ if(d02rho==1 && (matchK[3]>=0 && matchPi[0]>=0 && matchPi[1]>=0 && matchPi[2]>=0)) {*hyp2 = 1; output = 1;} //d02 hyp
+ if(d0bar1rho==1 && (matchK[0]>=0 && matchPi[1]>=0 && matchPi[2]>=0 && matchPi[3]>=0)) {*hyp3 = 1; output = 1;} //d0bar1 hyp
+ if(d0bar2rho==1 && (matchK[2]>=0 && matchPi[0]>=0 && matchPi[1]>=0 && matchPi[3]>=0)) {*hyp4 = 1; output = 1;} //d0bar2 hyp
+
+ return output;\r
}\r
//---------------------------------------------------------------------------\r
+Int_t AliRDHFCutsD0toKpipipi::D01Selected(TObject* obj,Int_t selectionLevel) {\r
+ //\r
+ // Apply selection\r
+ //\r
+\r
+ if(!fCutsRD){\r
+ cout<<"Cut matrix not inizialized. Exit..."<<endl;\r
+ return 0;\r
+ }\r
+ //PrintAll();\r
+ AliAODRecoDecayHF4Prong* d=(AliAODRecoDecayHF4Prong*)obj;\r
+\r
+ if(!d){\r
+ cout<<"AliAODRecoDecayHF4Prong null"<<endl;\r
+ return 0;\r
+ }\r
+\r
+ Int_t returnvalue=0;\r
+\r
+ // selection on candidate\r
+ if(selectionLevel==AliRDHFCuts::kAll || \r
+ selectionLevel==AliRDHFCuts::kCandidate) {\r
+\r
+ Int_t ptbin=PtBin(d->Pt());\r
+ \r
+ Double_t mD0[2];\r
+ Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
+\r
+ d->InvMassD0(mD0);\r
+ if(TMath::Abs(mD0[0]-mD0PDG) < fCutsRD[GetGlobalIndex(0,ptbin)]) returnvalue = 1;\r
+ }\r
+\r
+ return returnvalue;\r
+}
+//---------------------------------------------------------------------------\r
+Int_t AliRDHFCutsD0toKpipipi::D02Selected(TObject* obj,Int_t selectionLevel) {\r
+ //\r
+ // Apply selection\r
+ //\r
+\r
+ if(!fCutsRD){\r
+ cout<<"Cut matrix not inizialized. Exit..."<<endl;\r
+ return 0;\r
+ }\r
+ //PrintAll();\r
+ AliAODRecoDecayHF4Prong* d=(AliAODRecoDecayHF4Prong*)obj;\r
+\r
+ if(!d){\r
+ cout<<"AliAODRecoDecayHF4Prong null"<<endl;\r
+ return 0;\r
+ }\r
+\r
+ Int_t returnvalue=0;\r
+\r
+ // selection on candidate\r
+ if(selectionLevel==AliRDHFCuts::kAll || \r
+ selectionLevel==AliRDHFCuts::kCandidate) {\r
+\r
+ Int_t ptbin=PtBin(d->Pt());\r
+ \r
+ Double_t mD0[2];\r
+ Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
+\r
+ d->InvMassD0(mD0);\r
+ if(TMath::Abs(mD0[1]-mD0PDG) < fCutsRD[GetGlobalIndex(0,ptbin)]) returnvalue = 1;
+ }\r
+\r
+ return returnvalue;\r
+}\r//---------------------------------------------------------------------------\r
+Int_t AliRDHFCutsD0toKpipipi::D0bar1Selected(TObject* obj,Int_t selectionLevel) {\r
+ //\r
+ // Apply selection\r
+ //\r
+\r
+ if(!fCutsRD){\r
+ cout<<"Cut matrix not inizialized. Exit..."<<endl;\r
+ return 0;\r
+ }\r
+ //PrintAll();\r
+ AliAODRecoDecayHF4Prong* d=(AliAODRecoDecayHF4Prong*)obj;\r
+\r
+ if(!d){\r
+ cout<<"AliAODRecoDecayHF4Prong null"<<endl;\r
+ return 0;\r
+ }\r
+\r
+ Int_t returnvalue=0;\r
+\r
+ // selection on candidate\r
+ if(selectionLevel==AliRDHFCuts::kAll || \r
+ selectionLevel==AliRDHFCuts::kCandidate) {\r
+\r
+ Int_t ptbin=PtBin(d->Pt());\r
+ \r
+ Double_t mD0bar[2];\r
+ Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
+
+ d->InvMassD0bar(mD0bar);\r
+ if(TMath::Abs(mD0bar[0]-mD0PDG) < fCutsRD[GetGlobalIndex(0,ptbin)]) returnvalue = 1;\r
+ }\r
+\r
+ return returnvalue;\r
+}
+//---------------------------------------------------------------------------\r
+Int_t AliRDHFCutsD0toKpipipi::D0bar2Selected(TObject* obj,Int_t selectionLevel) {\r
+ //\r
+ // Apply selection\r
+ //\r
+\r
+ if(!fCutsRD){\r
+ cout<<"Cut matrix not inizialized. Exit..."<<endl;\r
+ return 0;\r
+ }\r
+ //PrintAll();\r
+ AliAODRecoDecayHF4Prong* d=(AliAODRecoDecayHF4Prong*)obj;\r
+\r
+ if(!d){\r
+ cout<<"AliAODRecoDecayHF4Prong null"<<endl;\r
+ return 0;\r
+ }\r
+\r
+ Int_t returnvalue=0;\r
+\r
+ // selection on candidate\r
+ if(selectionLevel==AliRDHFCuts::kAll || \r
+ selectionLevel==AliRDHFCuts::kCandidate) {\r
+\r
+ Int_t ptbin=PtBin(d->Pt());\r
+ \r
+ Double_t mD0bar[2];\r
+ Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
+
+ d->InvMassD0bar(mD0bar);\r
+ if(TMath::Abs(mD0bar[1]-mD0PDG) < fCutsRD[GetGlobalIndex(0,ptbin)]) returnvalue = 1;\r
+ }\r
+\r
+ return returnvalue;\r
+}
+//---------------------------------------------------------------------------
+Bool_t AliRDHFCutsD0toKpipipi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
+{
+ //
+ // Checking if D0 is in fiducial acceptance region
+ //
+
+ if(pt > 5.) {
+ // applying cut for pt > 5 GeV
+ AliDebug(4,Form("pt of D0 = %f (> 5), cutting at |y| < 0.8\n",pt));
+ if (TMath::Abs(y) > 0.8){
+ return kFALSE;
+ }
+ } else {
+ // appliying smooth cut for pt < 5 GeV
+ Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
+ Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
+ AliDebug(4,Form("pt of D0 = %f (< 5), cutting according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
+ if (y < minFiducialY || y > maxFiducialY){
+ return kFALSE;
+ }
+ }
+
+ return kTRUE;
+}
+
--- /dev/null
+//gSystem->AddIncludePath("-I/$ALICE_ROOT/PWG3/vertexingHF");
+
+#include <Riostream.h>
+#include <TFile.h>
+//#include <AliRDHFCutsD0toKpipipi.h>
+#include <TClonesArray.h>
+#include <TParameter.h>
+
+//Author: Chiara Bianchin, cbianchi@pd.infn.it, Fabio Colamaria, fabio.colamaria@ba.infn.it
+
+//macro to make a .root file which contains an AliRDHFCutsD0toKpipipi for AliAnalysisTaskSESelectHF4Prong task
+
+void MakeCuts4Charm4Prong(){
+
+gSystem->Load("libTree.so");
+gSystem->Load("libGeom.so");
+gSystem->Load("libPhysics.so");
+gSystem->Load("libVMC.so");
+gSystem->Load("libMinuit.so");
+gSystem->Load("libSTEERBase.so");
+gSystem->Load("libESD.so");
+gSystem->Load("libAOD.so");
+gSystem->Load("libANALYSIS.so");
+gSystem->Load("libANALYSISalice.so");
+gSystem->Load("libCORRFW.so");
+gSystem->Load("libPWG3base.so");
+gSystem->Load("libPWG3vertexingHF.so");
+gSystem->Load("libPWG3muon.so");
+
+ AliRDHFCutsD0toKpipipi* RDHFCharm4Prong=new AliRDHFCutsD0toKpipipi();
+ RDHFCharm4Prong->SetName("Charm4ProngCuts");
+ RDHFCharm4Prong->SetTitle("Cuts for D0 analysis");
+
+ AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
+ esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
+ //default
+ esdTrackCuts->SetRequireTPCRefit(kTRUE);
+ esdTrackCuts->SetRequireITSRefit(kTRUE);
+ esdTrackCuts->SetMinNClustersITS(4); // default is 5
+ //esdTrackCuts->SetMinNClustersTPC(70);
+ esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
+ AliESDtrackCuts::kAny);
+ // default is kBoth, otherwise kAny
+ esdTrackCuts->SetMinDCAToVertexXY(0.);
+ esdTrackCuts->SetPtRange(0.3,1.e10);
+
+ RDHFCharm4Prong->AddTrackCuts(esdTrackCuts);
+
+
+
+ const Int_t nvars=9;
+
+ const Int_t nptbins=5;
+ Float_t* ptbins;
+ ptbins=new Float_t[nptbins+1];
+ ptbins[0]=0.;
+ ptbins[1]=4.5;
+ ptbins[2]=6.;
+ ptbins[3]=8.;
+ ptbins[4]=12.;
+ ptbins[5]=16.;
+
+ RDHFCharm4Prong->SetPtBins(nptbins+1,ptbins);
+
+
+ Float_t** rdcutsvalmine;
+ rdcutsvalmine=new Float_t*[nvars];
+ for(Int_t iv=0;iv<nvars;iv++){
+ rdcutsvalmine[iv]=new Float_t[nptbins];
+ }
+
+ //setting my cut values
+ //cuts order
+ // printf(" |M-MD0| [GeV]");
+ // printf(" DCA [cm]");
+ // printf(" Dist 2-trk Vtx to PrimVtx [cm]");
+ // printf(" Dist 3-trk Vtx to PrimVtx [cm]");
+ // printf(" Dist 4-trk Vtx to PrimVtx [cm]");
+ // printf(" cosThetaPoint");
+ // printf(" pt [GeV/c]");
+ // printf(" rho mass");
+ // printf(" PID cut");
+
+ //setting my cut values
+ rdcutsvalmine[0][0]=0.2;
+ rdcutsvalmine[1][0]=0.025;
+ rdcutsvalmine[2][0]=0.0900;
+ rdcutsvalmine[3][0]=0.0900;
+ rdcutsvalmine[4][0]=0.0600;
+ rdcutsvalmine[5][0]=0.995;
+ rdcutsvalmine[6][0]=0.;
+ rdcutsvalmine[7][0]=0.1;
+ rdcutsvalmine[8][0]=0.;
+
+ rdcutsvalmine[0][1]=0.2;
+ rdcutsvalmine[1][1]=0.015;
+ rdcutsvalmine[2][1]=0.1200;
+ rdcutsvalmine[3][1]=0.1000;
+ rdcutsvalmine[4][1]=0.0850;
+ rdcutsvalmine[5][1]=0.99;
+ rdcutsvalmine[6][1]=0.;
+ rdcutsvalmine[7][1]=0.1;
+ rdcutsvalmine[8][1]=0.;
+
+ rdcutsvalmine[0][2]=0.2;
+ rdcutsvalmine[1][2]=0.025;
+ rdcutsvalmine[2][2]=0.0975;
+ rdcutsvalmine[3][2]=0.0900;
+ rdcutsvalmine[4][2]=0.0800;
+ rdcutsvalmine[5][2]=0.99;
+ rdcutsvalmine[6][2]=0.;
+ rdcutsvalmine[7][2]=0.1;
+ rdcutsvalmine[8][2]=0.;
+
+ rdcutsvalmine[0][3]=0.2;
+ rdcutsvalmine[1][3]=0.015;
+ rdcutsvalmine[2][3]=0.0975;
+ rdcutsvalmine[3][3]=0.0700;
+ rdcutsvalmine[4][3]=0.0750;
+ rdcutsvalmine[5][3]=0.9825;
+ rdcutsvalmine[6][3]=0.;
+ rdcutsvalmine[7][3]=0.1;
+ rdcutsvalmine[8][3]=0.;
+
+ rdcutsvalmine[0][4]=rdcutsvalmine[0][3];
+ rdcutsvalmine[1][4]=rdcutsvalmine[1][3];
+ rdcutsvalmine[2][4]=rdcutsvalmine[2][3];
+ rdcutsvalmine[3][4]=rdcutsvalmine[3][3];
+ rdcutsvalmine[4][4]=rdcutsvalmine[4][3];
+ rdcutsvalmine[5][4]=rdcutsvalmine[5][3];
+ rdcutsvalmine[6][4]=rdcutsvalmine[6][3];
+ rdcutsvalmine[7][4]=rdcutsvalmine[7][3];
+ rdcutsvalmine[8][4]=rdcutsvalmine[8][3];
+
+ RDHFCharm4Prong->SetCuts(nvars,nptbins,rdcutsvalmine);
+ cout<<"This is the odject I'm going to save:"<<endl;
+ RDHFCharm4Prong->PrintAll();
+ TFile* fout=new TFile("Charm4ProngCutsDef_16GeV.root","recreate"); //set this!!
+ fout->cd();
+ RDHFCharm4Prong->Write();
+ fout->Close();
+
+}
+
+void MakeCuts4Charm4ProngForMaxim(){
+
+gSystem->Clear();
+gSystem->Load("libTree.so");
+gSystem->Load("libGeom.so");
+gSystem->Load("libPhysics.so");
+gSystem->Load("libVMC.so");
+gSystem->Load("libMinuit.so");
+gSystem->Load("libSTEERBase.so");
+gSystem->Load("libESD.so");
+gSystem->Load("libAOD.so");
+gSystem->Load("libANALYSIS.so");
+gSystem->Load("libANALYSISalice.so");
+gSystem->Load("libCORRFW.so");
+gSystem->Load("libPWG3base.so");
+gSystem->Load("libPWG3vertexingHF.so");
+gSystem->Load("libPWG3muon.so");
+
+ AliRDHFCutsD0toKpipipi* RDHFCharm4Prong=new AliRDHFCutsD0toKpipipi();
+ RDHFCharm4Prong->SetName("loosercuts");
+ RDHFCharm4Prong->SetTitle("Cuts for significance maximization");
+
+ AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
+ esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
+ //default
+ esdTrackCuts->SetRequireTPCRefit(kTRUE);
+ esdTrackCuts->SetRequireITSRefit(kTRUE);
+ //esdTrackCuts->SetMinNClustersITS(4);
+ //esdTrackCuts->SetMinNClustersTPC(120);
+
+ esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
+ esdTrackCuts->SetMinDCAToVertexXY(0.);
+ esdTrackCuts->SetEtaRange(-0.8,0.8);
+ esdTrackCuts->SetPtRange(0.3,1.e10);
+
+ RDHFCharm4Prong->AddTrackCuts(esdTrackCuts);
+
+ const Int_t nvars=9;
+
+ const Int_t nptbins=5; //change this when adding pt bins!
+ Float_t ptbins[nptbins+1];
+ ptbins[0]=0.;
+ ptbins[1]=4.5;
+ ptbins[2]=6.;
+ ptbins[3]=8.;
+ ptbins[4]=12.;
+ ptbins[5]=16.; //o 25!
+
+ RDHFCharm4Prong->SetPtBins(nptbins+1,ptbins);
+
+ Float_t** rdcutsvalmine;
+ rdcutsvalmine=new Float_t*[nvars];
+ for(Int_t iv=0;iv<nvars;iv++){
+ rdcutsvalmine[iv]=new Float_t[nptbins];
+ }
+
+ //setting my cut values
+ //cuts order
+ // printf(" |M-MD0| [GeV]");
+ // printf(" DCA [cm]");
+ // printf(" Dist 2-trk Vtx to PrimVtx [cm]");
+ // printf(" Dist 3-trk Vtx to PrimVtx [cm]");
+ // printf(" Dist 4-trk Vtx to PrimVtx [cm]");
+ // printf(" cosThetaPoint");
+ // printf(" pt [GeV/c]");
+ // printf(" rho mass");
+ // printf(" PID cut");
+
+ Float_t cutsMatrixD0toKpipipiStand[nptbins][nvars]={{0.2,0.04,0.04,0.04,0.04,0.90,0.,0.1,0.},
+
+{0.2,0.04,0.04,0.04,0.04,0.90,0.,0.1,0.},
+
+{0.2,0.04,0.04,0.04,0.04,0.90,0.,0.1,0.},
+
+{0.2,0.04,0.04,0.04,0.04,0.90,0.,0.1,0.},
+
+{0.2,0.04,0.04,0.04,0.04,0.90,0.,0.1,0.}};
+
+ //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
+ Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
+ for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
+ for (Int_t ibin=0;ibin<nptbins;ibin++){
+ for (Int_t ivar = 0; ivar<nvars; ivar++){
+ cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpipipiStand[ibin][ivar];
+ }
+ }
+ RDHFCharm4Prong->SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
+
+ Int_t nvarsforopt=RDHFCharm4Prong->GetNVarsForOpt();
+ Int_t dim=5; //set this!!
+ Bool_t *boolforopt;
+ boolforopt=new Bool_t[nvars];
+ if(dim>nvarsforopt){
+ cout<<"Number of variables for optimization has probably changed, check and edit accordingly"<<endl;
+ return;
+ } else {
+ if(dim==nvarsforopt){
+ boolforopt=RDHFCharm4Prong->GetVarsForOpt();
+ }else{
+ TString *names;
+ names=new TString[nvars];
+ TString answer="";
+ Int_t checktrue=0;
+ names=RDHFCharm4Prong->GetVarNames();
+ for(Int_t i=0;i<nvars;i++){
+ cout<<names[i]<<" for opt? (y/n)"<<endl;
+ cin>>answer;
+ if(answer=="y") {
+ boolforopt[i]=kTRUE;
+ checktrue++;
+ }
+ else boolforopt[i]=kFALSE;
+ }
+ if (checktrue!=dim) {
+ cout<<"Error! You set "<<checktrue<<" kTRUE instead of "<<dim<<endl;
+ return;
+ }
+ RDHFCharm4Prong->SetVarsForOpt(dim,boolforopt);
+ }
+ }
+
+// Tight values for variables to be maximized
+// (DCA, Dist12, Dist3, Dist4, CosThetaPoint)
+// Indexes: variable, ptbin
+
+ Float_t tighterval[5][5];
+
+ //number of steps for each variable is set in the AddTask arguments (default=8)
+
+ tighterval[0][0]=0.025;
+ tighterval[1][0]=0.09;
+ tighterval[2][0]=0.09;
+ tighterval[3][0]=0.06;
+ tighterval[4][0]=0.955;
+
+ tighterval[0][1]=0.015;
+ tighterval[1][1]=0.12;
+ tighterval[2][1]=0.10;
+ tighterval[3][1]=0.085;
+ tighterval[4][1]=0.99;
+
+ tighterval[0][2]=0.025;
+ tighterval[1][2]=0.0975;
+ tighterval[2][2]=0.09;
+ tighterval[3][2]=0.08;
+ tighterval[4][2]=0.99;
+
+ tighterval[0][3]=0.015;
+ tighterval[1][3]=0.0975;
+ tighterval[2][3]=0.07;
+ tighterval[3][3]=0.075;
+ tighterval[4][3]=0.9825;
+
+ tighterval[0][4]=0.015;
+ tighterval[1][4]=0.0975;
+ tighterval[2][4]=0.07;
+ tighterval[3][4]=0.075;
+ tighterval[4][4]=0.9825;
+
+ TString name="";
+ Int_t arrdim=dim*nptbins;
+ cout<<"Will save "<<arrdim<<" TParameter<float>"<<endl;
+ TClonesArray max("TParameter<float>",arrdim);
+ for(Int_t ival=0;ival<dim;ival++){
+ for(Int_t jpt=0;jpt<nptbins;jpt++){
+ name=Form("par%dptbin%d",ival,jpt);
+ cout<<"Setting "<<name.Data()<<" to "<<tighterval[ival][jpt]<<endl;
+ new(max[jpt*dim+ival])TParameter<float>(name.Data(),tighterval[ival][jpt]);
+ }
+ }
+
+ Bool_t flagPID=kFALSE;
+ RDHFCharm4Prong->SetUsePID(flagPID);
+
+ RDHFCharm4Prong->PrintAll();
+ printf("Use PID? %s\n",flagPID ? "yes" : "no");
+
+/*
+ //pid settings
+ AliAODPidHF* pidObj=new AliAODPidHF();
+ //pidObj->SetName("pid4D0");
+ Int_t mode=1;
+ const Int_t nlims=2;
+ Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
+ Bool_t compat=kTRUE; //effective only for this mode
+ Bool_t asym=kTRUE;
+ Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
+ pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
+ pidObj->SetMatch(mode);
+ pidObj->SetPLimit(plims,nlims);
+ pidObj->SetSigma(sigmas);
+ pidObj->SetCompat(compat);
+ pidObj->SetTPC(kTRUE);
+ pidObj->SetTOF(kTRUE);
+ RDHFCharm4Prong->SetPidHF(pidObj);
+
+ RDHFCharm4Prong->SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
+*/
+
+ //activate pileup rejection (for pp)
+ RDHFCharm4Prong->SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
+
+ //Do not recalculate the vertex
+ RDHFCharm4Prong->SetRemoveDaughtersFromPrim(kTRUE); //activate for pp
+
+/* Per il pp levo la centralità!
+//Metto kCentOff! Per il Pb c'era kCentV0M
+ TString cent="";
+ //centrality selection (Pb-Pb)
+ Float_t minc=20,maxc=80;
+ RDHFCharm4Prong->SetMinCentrality(minc);
+ RDHFCharm4Prong->SetMaxCentrality(maxc);
+ cent=Form("%.0f%.0f",minc,maxc);
+ RDHFCharm4Prong->SetUseCentrality(AliRDHFCuts::kCentOff); //kCentOff,kCentV0M,kCentTRK,kCentTKL,kCentCL1,kCentInvalid
+*/
+
+ //temporary
+// RDHFCharm4Prong->SetFixRefs();
+
+ TFile* fout=new TFile("Charm4ProngCutsForMaxim.root","recreate");
+ fout->cd();
+ RDHFCharm4Prong->Write();
+ max.Write();
+ fout->Close();
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff