PWG3/vertexingHF/AliAnalysisTaskSESelectHF4Prong.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /////////////////////////////////////////////////////////////
  17 //
  18 // AliAnalysisTaskSE for the selection of heavy flavor
  19 // decay candidates and creation a stand-alone AOD for
  20 // 4prong D0 decay.
  21 //
  22 // Author: A.Dainese, andrea.dainese@lnl.infn.it
  23 //         F.Colamaria, fabio.colamaria@ba.infn.it
  24 /////////////////////////////////////////////////////////////
  25
  26
  27 #include "Riostream.h"
  28 #include "TFile.h"
  29 #include "TList.h"
  30 #include "TH1F.h"
  31 #include "TH2F.h"
  32 #include "TClonesArray.h"
  33 #include "TDatabasePDG.h"
  34 #include "TROOT.h"
  35 #include "TCanvas.h"
  36 #include "TBits.h"
  37 #include "TNtuple.h"
  38
  39 #include "AliAnalysisDataSlot.h"
  40 #include "AliAnalysisDataContainer.h"
  41 #include "AliAnalysisManager.h"
  42 #include "AliAODHandler.h"
  43 #include "AliAODEvent.h"
  44 #include "AliAODVertex.h"
  45 #include "AliAODTrack.h"
  46 #include "AliAODRecoDecayHF4Prong.h"
  47 #include "AliAnalysisVertexingHF.h"
  48 #include "AliAnalysisTaskSE.h"
  49 #include "AliAnalysisTaskSESelectHF4Prong.h"
  50 #include "AliAODPidHF.h"
  51 #include "AliRDHFCuts.h"
  52
  53 ClassImp(AliAnalysisTaskSESelectHF4Prong)
  54
  55
  56 //________________________________________________________________________
  57 AliAnalysisTaskSESelectHF4Prong::AliAnalysisTaskSESelectHF4Prong():
  58 AliAnalysisTaskSE(),
  59 fVerticesHFTClArr(0),
  60 fCharm4ProngTClArr(0),
  61 fOutput(0),
  62 fOutput2(0),
  63 fOutput3(0),
  64 fOutput4(0),
  65 fOutput5(0),
  66 fOutputC(0),
  67 fhInvMassD0Sum_10Mev_Bin1(0),
  68 fhInvMassD0barSum_10Mev_Bin1(0),
  69 fhInvMassSumAll_10Mev_Bin1(0),
  70 fhInvMassD0Sum_5Mev_Bin1(0),
  71 fhInvMassD0barSum_5Mev_Bin1(0),
  72 fhInvMassSumAll_5Mev_Bin1(0),
  73 fhInvMassD0Sum_10Mev_Bin2(0),
  74 fhInvMassD0barSum_10Mev_Bin2(0),
  75 fhInvMassSumAll_10Mev_Bin2(0),
  76 fhInvMassD0Sum_5Mev_Bin2(0),
  77 fhInvMassD0barSum_5Mev_Bin2(0),
  78 fhInvMassSumAll_5Mev_Bin2(0),
  79 fhInvMassD0Sum_10Mev_Bin3(0),
  80 fhInvMassD0barSum_10Mev_Bin3(0),
  81 fhInvMassSumAll_10Mev_Bin3(0),
  82 fhInvMassD0Sum_5Mev_Bin3(0),
  83 fhInvMassD0barSum_5Mev_Bin3(0),
  84 fhInvMassSumAll_5Mev_Bin3(0),
  85 fhInvMassD0Sum_10Mev_Bin4(0),
  86 fhInvMassD0barSum_10Mev_Bin4(0),
  87 fhInvMassSumAll_10Mev_Bin4(0),
  88 fhInvMassD0Sum_5Mev_Bin4(0),
  89 fhInvMassD0barSum_5Mev_Bin4(0),
  90 fhInvMassSumAll_5Mev_Bin4(0),
  91 fhInvMassD0Sum_10Mev_Bin5(0),
  92 fhInvMassD0barSum_10Mev_Bin5(0),
  93 fhInvMassSumAll_10Mev_Bin5(0),
  94 fhInvMassD0Sum_5Mev_Bin5(0),
  95 fhInvMassD0barSum_5Mev_Bin5(0),
  96 fhInvMassSumAll_5Mev_Bin5(0),
  97 fhInvMassMultipleOnly_Bin1(0),
  98 fhInvMassMultipleOnly_Bin2(0),
  99 fhInvMassMultipleOnly_Bin3(0),
 100 fhInvMassMultipleOnly_Bin4(0),
 101 fhInvMassMultipleOnly_Bin5(0),
 102 fScatterP4PID(0),
 103 fPtVsY(0),
 104 fPtVsYAll(0),
 105 fEventCounter(0),
 106 fCutDCA(0),
 107 fCutDCA3(0),
 108 fCutDCA2(0),
 109 fCutDCA5(0),
 110 fCutVertexDist2(0),
 111 fCutVertexDist3(0),
 112 fCutVertexDist4(0),
 113 fCutCosinePoint(0),
 114 fCutPt(0),
 115 fCutY(0),
 116 fPIDSel(0),
 117 fPIDSel_Bin1(0),
 118 fPIDSel_Bin2(0),
 119 fPIDSel_Bin3(0),
 120 fPIDSel_Bin4(0),
 121 fPIDSel_Bin5(0),
 122 fMultipleHyps(0),
 123 fMultipleHypsType(0),
 124 fPtSel(0),
 125 fCuts(0),
 126 fVHF(0)
 127 {
 128   // Default constructor
 129 }
 130
 131 //________________________________________________________________________
 132 AliAnalysisTaskSESelectHF4Prong::AliAnalysisTaskSESelectHF4Prong(const char *name,AliRDHFCutsD0toKpipipi* cuts):
 133 AliAnalysisTaskSE(name),
 134 fVerticesHFTClArr(0),
 135 fCharm4ProngTClArr(0),
 136 fOutput(0),
 137 fOutput2(0),
 138 fOutput3(0),
 139 fOutput4(0),
 140 fOutput5(0),
 141 fOutputC(0),
 142 fhInvMassD0Sum_10Mev_Bin1(0),
 143 fhInvMassD0barSum_10Mev_Bin1(0),
 144 fhInvMassSumAll_10Mev_Bin1(0),
 145 fhInvMassD0Sum_5Mev_Bin1(0),
 146 fhInvMassD0barSum_5Mev_Bin1(0),
 147 fhInvMassSumAll_5Mev_Bin1(0),
 148 fhInvMassD0Sum_10Mev_Bin2(0),
 149 fhInvMassD0barSum_10Mev_Bin2(0),
 150 fhInvMassSumAll_10Mev_Bin2(0),
 151 fhInvMassD0Sum_5Mev_Bin2(0),
 152 fhInvMassD0barSum_5Mev_Bin2(0),
 153 fhInvMassSumAll_5Mev_Bin2(0),
 154 fhInvMassD0Sum_10Mev_Bin3(0),
 155 fhInvMassD0barSum_10Mev_Bin3(0),
 156 fhInvMassSumAll_10Mev_Bin3(0),
 157 fhInvMassD0Sum_5Mev_Bin3(0),
 158 fhInvMassD0barSum_5Mev_Bin3(0),
 159 fhInvMassSumAll_5Mev_Bin3(0),
 160 fhInvMassD0Sum_10Mev_Bin4(0),
 161 fhInvMassD0barSum_10Mev_Bin4(0),
 162 fhInvMassSumAll_10Mev_Bin4(0),
 163 fhInvMassD0Sum_5Mev_Bin4(0),
 164 fhInvMassD0barSum_5Mev_Bin4(0),
 165 fhInvMassSumAll_5Mev_Bin4(0),
 166 fhInvMassD0Sum_10Mev_Bin5(0),
 167 fhInvMassD0barSum_10Mev_Bin5(0),
 168 fhInvMassSumAll_10Mev_Bin5(0),
 169 fhInvMassD0Sum_5Mev_Bin5(0),
 170 fhInvMassD0barSum_5Mev_Bin5(0),
 171 fhInvMassSumAll_5Mev_Bin5(0),
 172 fhInvMassMultipleOnly_Bin1(0),
 173 fhInvMassMultipleOnly_Bin2(0),
 174 fhInvMassMultipleOnly_Bin3(0),
 175 fhInvMassMultipleOnly_Bin4(0),
 176 fhInvMassMultipleOnly_Bin5(0),
 177 fScatterP4PID(0),
 178 fPtVsY(0),
 179 fPtVsYAll(0),
 180 fEventCounter(0),
 181 fCutDCA(0),
 182 fCutDCA3(0),
 183 fCutDCA2(0),
 184 fCutDCA5(0),
 185 fCutVertexDist2(0),
 186 fCutVertexDist3(0),
 187 fCutVertexDist4(0),
 188 fCutCosinePoint(0),
 189 fCutPt(0),
 190 fCutY(0),
 191 fPIDSel(0),
 192 fPIDSel_Bin1(0),
 193 fPIDSel_Bin2(0),
 194 fPIDSel_Bin3(0),
 195 fPIDSel_Bin4(0),
 196 fPIDSel_Bin5(0),
 197 fMultipleHyps(0),
 198 fMultipleHypsType(0),
 199 fPtSel(0),
 200 fCuts(0),
 201 fVHF(0)
 202 {
 203   // Standard constructor
 204
 205   fCuts=cuts;
 206
 207   // Input slot #0 works with an Ntuple
 208 //  DefineInput(0, TTree::Class());
 209
 210   // Output slot #0 writes into a TTree container
 211   // Output slots #1-6 writes into a TList container
 212   DefineOutput(0, TTree::Class()); //default
 213   DefineOutput(1, TList::Class()); //histos inv. mass bin1
 214   DefineOutput(2, TList::Class()); //histos inv. mass bin2
 215   DefineOutput(3, TList::Class()); //histos inv. mass bin3
 216   DefineOutput(4, TList::Class()); //histos inv. mass bin4
 217   DefineOutput(5, TList::Class()); //histos inv. mass bin5
 218   DefineOutput(6, TList::Class()); //histos of cuts
 219   DefineOutput(7, AliRDHFCutsD0toKpipipi::Class()); //cuts
 220 }
 221
 222 //________________________________________________________________________
 223 AliAnalysisTaskSESelectHF4Prong::~AliAnalysisTaskSESelectHF4Prong()
 224 {
 225   // Destructor
 226
 227   if (fOutput) {
 228     delete fOutput;
 229     fOutput = 0;
 230   }
 231   if (fOutput2) {
 232     delete fOutput2;
 233     fOutput2 = 0;
 234   }
 235   if (fOutput3) {
 236     delete fOutput3;
 237     fOutput3 = 0;
 238   }
 239   if (fOutput4) {
 240     delete fOutput4;
 241     fOutput4 = 0;
 242   }
 243   if (fOutput5) {
 244     delete fOutput5;
 245     fOutput5 = 0;
 246   }
 247   if (fOutputC) {
 248     delete fOutputC;
 249     fOutputC = 0;
 250   }
 251   if (fCuts) {
 252     delete fCuts;
 253     fCuts = 0;
 254   }
 255   if (fVHF) {
 256     delete fVHF;
 257     fVHF = 0;
 258   }
 259
 260 }
 261
 262 //________________________________________________________________________
 263 void AliAnalysisTaskSESelectHF4Prong::Init()
 264 {
 265   // Initialization
 266
 267   if(fDebug > 1) printf("AnalysisTaskSESelectHF4Prong::Init() \n");
 268
 269   return;
 270 }
 271
 272 //________________________________________________________________________
 273 void AliAnalysisTaskSESelectHF4Prong::UserCreateOutputObjects()
 274 {
 275   // Create the output container
 276   //
 277   if(fDebug > 1) printf("AnalysisTaskSESelectHF4Prong::UserCreateOutputObjects() \n");
 278
 279   fVerticesHFTClArr = new TClonesArray("AliAODVertex", 0);
 280   fVerticesHFTClArr->SetName("VerticesHF");
 281   AddAODBranch("TClonesArray", &fVerticesHFTClArr);
 282
 283   fCharm4ProngTClArr = new TClonesArray("AliAODRecoDecayHF4Prong", 0);
 284   fCharm4ProngTClArr->SetName("Charm4Prong");
 285   AddAODBranch("TClonesArray", &fCharm4ProngTClArr);
 286
 287   fOutput = new TList();
 288   fOutput->SetOwner();
 289
 290   fOutput2 = new TList();
 291   fOutput2->SetOwner();
 292
 293   fOutput3 = new TList();
 294   fOutput3->SetOwner();
 295
 296   fOutput4 = new TList();
 297   fOutput4->SetOwner();
 298
 299   fOutput5 = new TList();
 300   fOutput5->SetOwner();
 301
 302   fOutputC = new TList();
 303   fOutputC->SetOwner();
 304
 305   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);
 306   fhInvMassD0Sum_10Mev_Bin1->Sumw2(); //Create structure to store sum of squares of weights
 307   fhInvMassD0Sum_10Mev_Bin1->SetMinimum(0);
 308   fOutput->Add(fhInvMassD0Sum_10Mev_Bin1);
 309
 310   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);
 311   fhInvMassD0barSum_10Mev_Bin1->Sumw2(); //Create structure to store sum of squares of weights
 312   fhInvMassD0barSum_10Mev_Bin1->SetMinimum(0);
 313   fOutput->Add(fhInvMassD0barSum_10Mev_Bin1);
 314
 315   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);
 316   fhInvMassSumAll_10Mev_Bin1->Sumw2(); //Create structure to store sum of squares of weights
 317   fhInvMassSumAll_10Mev_Bin1->SetMinimum(0);
 318   fOutput->Add(fhInvMassSumAll_10Mev_Bin1);
 319
 320   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);
 321   fhInvMassD0Sum_5Mev_Bin1->Sumw2(); //Create structure to store sum of squares of weights
 322   fhInvMassD0Sum_5Mev_Bin1->SetMinimum(0);
 323   fOutput->Add(fhInvMassD0Sum_5Mev_Bin1);
 324
 325   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);
 326   fhInvMassD0barSum_5Mev_Bin1->Sumw2(); //Create structure to store sum of squares of weights
 327   fhInvMassD0barSum_5Mev_Bin1->SetMinimum(0);
 328   fOutput->Add(fhInvMassD0barSum_5Mev_Bin1);
 329
 330   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);
 331   fhInvMassSumAll_5Mev_Bin1->Sumw2(); //Create structure to store sum of squares of weights
 332   fhInvMassSumAll_5Mev_Bin1->SetMinimum(0);
 333   fOutput->Add(fhInvMassSumAll_5Mev_Bin1);
 334
 335   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);
 336   fhInvMassD0Sum_10Mev_Bin2->Sumw2(); //Create structure to store sum of squares of weights
 337   fhInvMassD0Sum_10Mev_Bin2->SetMinimum(0);
 338   fOutput2->Add(fhInvMassD0Sum_10Mev_Bin2);
 339
 340   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);
 341   fhInvMassD0barSum_10Mev_Bin2->Sumw2(); //Create structure to store sum of squares of weights
 342   fhInvMassD0barSum_10Mev_Bin2->SetMinimum(0);
 343   fOutput2->Add(fhInvMassD0barSum_10Mev_Bin2);
 344
 345   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);
 346   fhInvMassSumAll_10Mev_Bin2->Sumw2(); //Create structure to store sum of squares of weights
 347   fhInvMassSumAll_10Mev_Bin2->SetMinimum(0);
 348   fOutput2->Add(fhInvMassSumAll_10Mev_Bin2);
 349
 350   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);
 351   fhInvMassD0Sum_5Mev_Bin2->Sumw2(); //Create structure to store sum of squares of weights
 352   fhInvMassD0Sum_5Mev_Bin2->SetMinimum(0);
 353   fOutput2->Add(fhInvMassD0Sum_5Mev_Bin2);
 354
 355   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);
 356   fhInvMassD0barSum_5Mev_Bin2->Sumw2(); //Create structure to store sum of squares of weights
 357   fhInvMassD0barSum_5Mev_Bin2->SetMinimum(0);
 358   fOutput2->Add(fhInvMassD0barSum_5Mev_Bin2);
 359
 360   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);
 361   fhInvMassSumAll_5Mev_Bin2->Sumw2(); //Create structure to store sum of squares of weights
 362   fhInvMassSumAll_5Mev_Bin2->SetMinimum(0);
 363   fOutput2->Add(fhInvMassSumAll_5Mev_Bin2);
 364
 365   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);
 366   fhInvMassD0Sum_10Mev_Bin3->Sumw2(); //Create structure to store sum of squares of weights
 367   fhInvMassD0Sum_10Mev_Bin3->SetMinimum(0);
 368   fOutput3->Add(fhInvMassD0Sum_10Mev_Bin3);
 369
 370   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);
 371   fhInvMassD0barSum_10Mev_Bin3->Sumw2(); //Create structure to store sum of squares of weights
 372   fhInvMassD0barSum_10Mev_Bin3->SetMinimum(0);
 373   fOutput3->Add(fhInvMassD0barSum_10Mev_Bin3);
 374
 375   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);
 376   fhInvMassSumAll_10Mev_Bin3->Sumw2(); //Create structure to store sum of squares of weights
 377   fhInvMassSumAll_10Mev_Bin3->SetMinimum(0);
 378   fOutput3->Add(fhInvMassSumAll_10Mev_Bin3);
 379
 380   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);
 381   fhInvMassD0Sum_5Mev_Bin3->Sumw2(); //Create structure to store sum of squares of weights
 382   fhInvMassD0Sum_5Mev_Bin3->SetMinimum(0);
 383   fOutput3->Add(fhInvMassD0Sum_5Mev_Bin3);
 384
 385   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);
 386   fhInvMassD0barSum_5Mev_Bin3->Sumw2(); //Create structure to store sum of squares of weights
 387   fhInvMassD0barSum_5Mev_Bin3->SetMinimum(0);
 388   fOutput3->Add(fhInvMassD0barSum_5Mev_Bin3);
 389
 390   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);
 391   fhInvMassSumAll_5Mev_Bin3->Sumw2(); //Create structure to store sum of squares of weights
 392   fhInvMassSumAll_5Mev_Bin3->SetMinimum(0);
 393   fOutput3->Add(fhInvMassSumAll_5Mev_Bin3);
 394
 395   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);
 396   fhInvMassD0Sum_10Mev_Bin4->Sumw2(); //Create structure to store sum of squares of weights
 397   fhInvMassD0Sum_10Mev_Bin4->SetMinimum(0);
 398   fOutput4->Add(fhInvMassD0Sum_10Mev_Bin4);
 399
 400   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);
 401   fhInvMassD0barSum_10Mev_Bin4->Sumw2(); //Create structure to store sum of squares of weights
 402   fhInvMassD0barSum_10Mev_Bin4->SetMinimum(0);
 403   fOutput4->Add(fhInvMassD0barSum_10Mev_Bin4);
 404
 405   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);
 406   fhInvMassSumAll_10Mev_Bin4->Sumw2(); //Create structure to store sum of squares of weights
 407   fhInvMassSumAll_10Mev_Bin4->SetMinimum(0);
 408   fOutput4->Add(fhInvMassSumAll_10Mev_Bin4);
 409
 410   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);
 411   fhInvMassD0Sum_5Mev_Bin4->Sumw2(); //Create structure to store sum of squares of weights
 412   fhInvMassD0Sum_5Mev_Bin4->SetMinimum(0);
 413   fOutput4->Add(fhInvMassD0Sum_5Mev_Bin4);
 414
 415   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);
 416   fhInvMassD0barSum_5Mev_Bin4->Sumw2(); //Create structure to store sum of squares of weights
 417   fhInvMassD0barSum_5Mev_Bin4->SetMinimum(0);
 418   fOutput4->Add(fhInvMassD0barSum_5Mev_Bin4);
 419
 420   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);
 421   fhInvMassSumAll_5Mev_Bin4->Sumw2(); //Create structure to store sum of squares of weights
 422   fhInvMassSumAll_5Mev_Bin4->SetMinimum(0);
 423   fOutput4->Add(fhInvMassSumAll_5Mev_Bin4);
 424
 425   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);
 426   fhInvMassD0Sum_10Mev_Bin5->Sumw2(); //Create structure to store sum of squares of weights
 427   fhInvMassD0Sum_10Mev_Bin5->SetMinimum(0);
 428   fOutput5->Add(fhInvMassD0Sum_10Mev_Bin5);
 429
 430   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);
 431   fhInvMassD0barSum_10Mev_Bin5->Sumw2(); //Create structure to store sum of squares of weights
 432   fhInvMassD0barSum_10Mev_Bin5->SetMinimum(0);
 433   fOutput5->Add(fhInvMassD0barSum_10Mev_Bin5);
 434
 435   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);
 436   fhInvMassSumAll_10Mev_Bin5->Sumw2(); //Create structure to store sum of squares of weights
 437   fhInvMassSumAll_10Mev_Bin5->SetMinimum(0);
 438   fOutput5->Add(fhInvMassSumAll_10Mev_Bin5);
 439
 440   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);
 441   fhInvMassD0Sum_5Mev_Bin5->Sumw2(); //Create structure to store sum of squares of weights
 442   fhInvMassD0Sum_5Mev_Bin5->SetMinimum(0);
 443   fOutput5->Add(fhInvMassD0Sum_5Mev_Bin5);
 444
 445   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);
 446   fhInvMassD0barSum_5Mev_Bin5->Sumw2(); //Create structure to store sum of squares of weights
 447   fhInvMassD0barSum_5Mev_Bin5->SetMinimum(0);
 448   fOutput5->Add(fhInvMassD0barSum_5Mev_Bin5);
 449
 450   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);
 451   fhInvMassSumAll_5Mev_Bin5->Sumw2(); //Create structure to store sum of squares of weights
 452   fhInvMassSumAll_5Mev_Bin5->SetMinimum(0);
 453   fOutput5->Add(fhInvMassSumAll_5Mev_Bin5);
 454
 455   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);
 456   fhInvMassMultipleOnly_Bin1->Sumw2(); //Create structure to store sum of squares of weights
 457   fhInvMassMultipleOnly_Bin1->SetMinimum(0);
 458   fOutput->Add(fhInvMassMultipleOnly_Bin1);
 459
 460   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);
 461   fhInvMassMultipleOnly_Bin2->Sumw2(); //Create structure to store sum of squares of weights
 462   fhInvMassMultipleOnly_Bin2->SetMinimum(0);
 463   fOutput2->Add(fhInvMassMultipleOnly_Bin2);
 464
 465   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);
 466   fhInvMassMultipleOnly_Bin3->Sumw2(); //Create structure to store sum of squares of weights
 467   fhInvMassMultipleOnly_Bin3->SetMinimum(0);
 468   fOutput3->Add(fhInvMassMultipleOnly_Bin3);
 469
 470   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);
 471   fhInvMassMultipleOnly_Bin4->Sumw2(); //Create structure to store sum of squares of weights
 472   fhInvMassMultipleOnly_Bin4->SetMinimum(0);
 473   fOutput4->Add(fhInvMassMultipleOnly_Bin4);
 474
 475   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);
 476   fhInvMassMultipleOnly_Bin5->Sumw2(); //Create structure to store sum of squares of weights
 477   fhInvMassMultipleOnly_Bin5->SetMinimum(0);
 478   fOutput5->Add(fhInvMassMultipleOnly_Bin5);
 479
 480   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.);
 481   fScatterP4PID->SetMinimum(0);
 482   fOutput->Add(fScatterP4PID);
 483
 484   fPtVsY = new TH2F("fPtVsY", "Pt vs Y PPR Sel. Candidates; Pt [GeV/c]; Y",250,0.,25.,300,-3.,3.);
 485   fPtVsY->SetMinimum(0);
 486   fOutputC->Add(fPtVsY);
 487
 488   fPtVsYAll = new TH2F("fPtVsYAll", "Pt vs Y All Candidates; Pt [GeV/c]; Y",250,0.,25.,300,-3.,3.);
 489   fPtVsYAll->SetMinimum(0);
 490   fOutputC->Add(fPtVsYAll);
 491
 492   fEventCounter = new TH1F("fEventCounter", "N° of total events; NA; Events",1,0.,1.);
 493   fEventCounter->SetMinimum(0);
 494   fOutputC->Add(fEventCounter);
 495
 496   fCutDCA = new TH1F("fCutDCA", "DCA of candidate (couple); DCA [cm]; Entries/micron",500,0.,0.05);
 497   fCutDCA->SetMinimum(0);
 498   fOutputC->Add(fCutDCA);
 499
 500   fCutDCA3 = new TH1F("fCutDCA3", "DCA of candidate (trips); DCA [cm]; Entries/micron",500,0.,0.05);
 501   fCutDCA3->SetMinimum(0);
 502   fOutputC->Add(fCutDCA3);
 503
 504   fCutDCA2 = new TH1F("fCutDCA2", "DCA of candidate (quads1); DCA [cm]; Entries/micron",500,0.,0.05);
 505   fCutDCA2->SetMinimum(0);
 506   fOutputC->Add(fCutDCA2);
 507
 508   fCutDCA5 = new TH1F("fCutDCA5", "DCA of candidate (quads2); DCA [cm]; Entries/micron",500,0.,0.05);
 509   fCutDCA5->SetMinimum(0);
 510   fOutputC->Add(fCutDCA5);
 511
 512   fCutVertexDist2 = new TH1F("fCutVertexDist2", "Distance Vtx doubl.-Primary Vtx; Distance [cm]; Entries/15 micron",500,0.,0.75);
 513   fCutVertexDist2->SetMinimum(0);
 514   fOutputC->Add(fCutVertexDist2);
 515
 516   fCutVertexDist3 = new TH1F("fCutVertexDist3", "Distance Vtx trips-Primary Vtx; Distance [cm]; Entries/10 micron",500,0.,0.5);
 517   fCutVertexDist3->SetMinimum(0);
 518   fOutputC->Add(fCutVertexDist3);
 519
 520   fCutVertexDist4 = new TH1F("fCutVertexDist4", "Distance Vtx quads-Primary Vtx; Distance [cm]; Entries/5 micron",500,0.,0.25);
 521   fCutVertexDist4->SetMinimum(0);
 522   fOutputC->Add(fCutVertexDist4);
 523
 524   fCutCosinePoint = new TH1F("fCutCosinePoint", "Cosine of angle of pointing; Cos(Thetapt.); Entries/10^(-3)",250,0.75,1.);
 525   fCutCosinePoint->SetMinimum(0);
 526   fOutputC->Add(fCutCosinePoint);
 527
 528   fCutPt = new TH1F("fCutPt", "Pt of candidate D0; Pt [GeV/c]; Entries/5 MeV",3000,0.,15.);
 529   fCutPt->SetMinimum(0);
 530   fOutputC->Add(fCutPt);
 531
 532   fCutY = new TH1F("fCutY", "Y of candidate D0; Pt [GeV/c]; Entries/5 MeV",900,-9.,9.);
 533   fCutY->SetMinimum(0);
 534   fOutputC->Add(fCutY);
 535
 536   fPIDSel = new TH1F("fPIDSel", "Ratio of D0 selected by PID for Correction",3,0.,3.);
 537   fPIDSel->SetMinimum(0);
 538   fPIDSel->GetXaxis()->SetBinLabel(1,"D0allhyp All");
 539   fPIDSel->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
 540   fPIDSel->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
 541
 542   fPIDSel_Bin1 = new TH1F("fPIDSel_Bin1", "Ratio of D0 selected by PID for Correction",3,0.,3.);
 543   fPIDSel_Bin1->SetMinimum(0);
 544   fPIDSel_Bin1->GetXaxis()->SetBinLabel(1,"D0allhyp All");
 545   fPIDSel_Bin1->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
 546   fPIDSel_Bin1->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
 547
 548   fPIDSel_Bin2 = new TH1F("fPIDSel_Bin1", "Ratio of D0 selected by PID for Correction",3,0.,3.);
 549   fPIDSel_Bin2->SetMinimum(0);
 550   fPIDSel_Bin2->GetXaxis()->SetBinLabel(1,"D0allhyp All");
 551   fPIDSel_Bin2->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
 552   fPIDSel_Bin2->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
 553
 554   fPIDSel_Bin3 = new TH1F("fPIDSel_Bin1", "Ratio of D0 selected by PID for Correction",3,0.,3.);
 555   fPIDSel_Bin3->SetMinimum(0);
 556   fPIDSel_Bin3->GetXaxis()->SetBinLabel(1,"D0allhyp All");
 557   fPIDSel_Bin3->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
 558   fPIDSel_Bin3->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
 559
 560   fPIDSel_Bin4 = new TH1F("fPIDSel_Bin1", "Ratio of D0 selected by PID for Correction",3,0.,3.);
 561   fPIDSel_Bin4->SetMinimum(0);
 562   fPIDSel_Bin4->GetXaxis()->SetBinLabel(1,"D0allhyp All");
 563   fPIDSel_Bin4->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
 564   fPIDSel_Bin4->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
 565
 566   fPIDSel_Bin5 = new TH1F("fPIDSel_Bin1", "Ratio of D0 selected by PID for Correction",3,0.,3.);
 567   fPIDSel_Bin5->SetMinimum(0);
 568   fPIDSel_Bin5->GetXaxis()->SetBinLabel(1,"D0allhyp All");
 569   fPIDSel_Bin5->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
 570   fPIDSel_Bin5->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
 571
 572   fMultipleHyps = new TH1F("fMultipleHyps", "N. of hyp. accepted for each candidate (accounted N. times)",8,0.,8.);
 573   fMultipleHyps->SetMinimum(0);
 574   fMultipleHyps->GetXaxis()->SetBinLabel(1,"1 (PPR)");
 575   fMultipleHyps->GetXaxis()->SetBinLabel(2,"2 (PPR)");
 576   fMultipleHyps->GetXaxis()->SetBinLabel(3,"3 (PPR)");
 577   fMultipleHyps->GetXaxis()->SetBinLabel(4,"4 (PPR)");
 578   fMultipleHyps->GetXaxis()->SetBinLabel(5,"1 (PPR+PID)");
 579   fMultipleHyps->GetXaxis()->SetBinLabel(6,"2 (PPR+PID)");
 580   fMultipleHyps->GetXaxis()->SetBinLabel(7,"3 (PPR+PID)");
 581   fMultipleHyps->GetXaxis()->SetBinLabel(8,"4 (PPR+PID)");
 582
 583   fMultipleHypsType = new TH1F("fMultipleHypsType", "Type of hyp. accepted for each candidate",8,0.,8.);
 584   fMultipleHypsType->SetMinimum(0);
 585   fMultipleHypsType->GetXaxis()->SetBinLabel(1,"D0");
 586   fMultipleHypsType->GetXaxis()->SetBinLabel(2,"D0bar");
 587   fMultipleHypsType->GetXaxis()->SetBinLabel(3,"2D0");
 588   fMultipleHypsType->GetXaxis()->SetBinLabel(4,"2D0bar");
 589   fMultipleHypsType->GetXaxis()->SetBinLabel(5,"D0+D0bar");
 590   fMultipleHypsType->GetXaxis()->SetBinLabel(6,"2D0+D0bar");
 591   fMultipleHypsType->GetXaxis()->SetBinLabel(7,"D0+2D0bar");
 592   fMultipleHypsType->GetXaxis()->SetBinLabel(8,"2D0+2D0bar");
 593
 594   fOutputC->Add(fMultipleHyps);
 595   fOutputC->Add(fMultipleHypsType);
 596
 597   fOutputC->Add(fPIDSel);
 598   fOutput->Add(fPIDSel_Bin1);
 599   fOutput2->Add(fPIDSel_Bin2);
 600   fOutput3->Add(fPIDSel_Bin3);
 601   fOutput4->Add(fPIDSel_Bin4);
 602   fOutput5->Add(fPIDSel_Bin5);
 603
 604   fPtSel = new TH1F("fPtSel", "Pt of candidates accepted; Pt [GeV/c]; Entries/10 MeV",2000,0.,20.);
 605   fPtSel->SetMinimum(0);
 606   fOutputC->Add(fPtSel);
 607
 608   return;
 609 }
 610
 611 //________________________________________________________________________
 612 void AliAnalysisTaskSESelectHF4Prong::UserExec(Option_t */*option*/)
 613 {
 614   // Execute analysis for current event:
 615   // heavy flavor candidates selection and histograms
 616
 617   AliAODEvent *aodIn = dynamic_cast<AliAODEvent*> (InputEvent());
 618
 619   TClonesArray *inputArrayCharm4Prong = 0;
 620
 621   if(!aodIn && AODEvent() && IsStandardAOD()) {
 622     // In case there is an AOD handler writing a standard AOD, use the AOD
 623     // event in memory rather than the input (ESD) event.
 624     aodIn = dynamic_cast<AliAODEvent*> (AODEvent());
 625     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 626     // have to taken from the AOD event hold by the AliAODExtension
 627     AliAODHandler* aodHandler = (AliAODHandler*)
 628       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 629     if(aodHandler->GetExtensions()) {
 630       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 631       AliAODEvent *aodFromExt = ext->GetAOD();
 632       // load D0 candidates
 633       inputArrayCharm4Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
 634     }
 635   } else {
 636     // load D0 candidates
 637     inputArrayCharm4Prong=(TClonesArray*)aodIn->GetList()->FindObject("Charm4Prong");
 638   }
 639
 640   if(!inputArrayCharm4Prong) {
 641     printf("AliAnalysisTaskSESelectHF4Prong::UserExec: D0to3Kpi branch not found!\n");
 642     return;
 643   }
 644
 645   //print event info
 646 //  aodIn->GetHeader()->Print();
 647
 648   //Event counter ++
 649   fEventCounter->Fill(0);
 650
 651   // fix for temporary bug in ESDfilter
 652   // the AODs with null vertex pointer didn't pass the PhysSel
 653   if(!aodIn->GetPrimaryVertex()) return;
 654
 655   // primary vertex
 656   AliAODVertex *vtx1 = (AliAODVertex*)aodIn->GetPrimaryVertex();
 657 //  vtx1->Print();
 658
 659   // make trkIDtoEntry register (temporary)
 660   Int_t trkIDtoEntry[100000];
 661   for(Int_t it=0;it<aodIn->GetNumberOfTracks();it++) {
 662     AliAODTrack *track = aodIn->GetTrack(it);
 663     trkIDtoEntry[track->GetID()]=it;
 664   }
 665
 666   Int_t iOutVerticesHF=0,iOutCharm4Prong=0;
 667   fVerticesHFTClArr->Delete();
 668   iOutVerticesHF = fVerticesHFTClArr->GetEntriesFast();
 669   TClonesArray &verticesHFRef = *fVerticesHFTClArr;
 670   fCharm4ProngTClArr->Delete();
 671   iOutCharm4Prong = fCharm4ProngTClArr->GetEntriesFast();
 672   TClonesArray &aodCharm4ProngRef = *fCharm4ProngTClArr;
 673
 674   // loop over D0->K3pi candidates
 675   Int_t nInCharm4Prong = inputArrayCharm4Prong->GetEntriesFast();
 676   printf("Number of D0->K3pi: %d\n",nInCharm4Prong);
 677
 678   for (Int_t iCharm4Prong = 0; iCharm4Prong < nInCharm4Prong; iCharm4Prong++) {
 679     AliAODRecoDecayHF4Prong *dIn = (AliAODRecoDecayHF4Prong*)inputArrayCharm4Prong->UncheckedAt(iCharm4Prong);
 680     Bool_t unsetvtx=kFALSE;
 681
 682     if(!dIn->GetOwnPrimaryVtx()) {
 683       dIn->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 684       unsetvtx=kTRUE;
 685     }
 686
 687    //fill histos of cuts
 688    Double_t dca = dIn->GetDCA();
 689    Double_t dist2 = dIn->GetDist12toPrim();
 690    Double_t dist3 = dIn->GetDist3toPrim();
 691    Double_t dist4 = dIn->GetDist4toPrim();
 692    Double_t cosine = dIn->CosPointingAngle();
 693    Double_t ptPart = dIn->Pt();
 694    Double_t yPart = dIn->YD0();
 695    Double_t dcatrip = dIn->GetDCA(3);
 696    Double_t dcaquad1 = dIn->GetDCA(2);
 697    Double_t dcaquad2 = dIn->GetDCA(5);
 698
 699    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
 700
 701    fCutDCA->Fill(dca);
 702    fCutDCA3->Fill(dcatrip);
 703    fCutDCA2->Fill(dcaquad1);
 704    fCutDCA5->Fill(dcaquad2);
 705    fCutVertexDist2->Fill(dist2);
 706    fCutVertexDist3->Fill(dist3);
 707    fCutVertexDist4->Fill(dist4);
 708    fCutCosinePoint->Fill(cosine);
 709    fCutPt->Fill(ptPart);
 710    fCutY->Fill(yPart);
 711    fPtVsYAll->Fill(ptPart,yPart);
 712
 713    //flags initialization
 714    fSelected = fCuts->IsSelected(dIn,AliRDHFCuts::kCandidate);
 715    Int_t selD01 = fCuts->D01Selected(dIn,AliRDHFCuts::kCandidate);
 716    Int_t selD02 = fCuts->D02Selected(dIn,AliRDHFCuts::kCandidate);
 717    Int_t selD0bar1 = fCuts->D0bar1Selected(dIn,AliRDHFCuts::kCandidate);
 718    Int_t selD0bar2 = fCuts->D0bar2Selected(dIn,AliRDHFCuts::kCandidate);
 719    Int_t flagAccLim = 1;
 720
 721    //Limited Acceptance
 722    if(ptPart > 5.) {
 723      if (TMath::Abs(yPart) > 0.8) flagAccLim = 0;
 724    }
 725    else {
 726      Double_t maxFiducialY = -0.2/15*ptPart*ptPart+1.9/15*ptPart+0.5;
 727      Double_t minFiducialY = 0.2/15*ptPart*ptPart-1.9/15*ptPart-0.5;
 728      if (yPart < minFiducialY || yPart > maxFiducialY) flagAccLim = 0;;
 729    }
 730
 731    //number of CANDIDATES (regardless of hypotheses) passing PPR
 732    if(fSelected==1||fSelected==2||fSelected==3) {
 733       fPtVsY->Fill(ptPart,yPart);
 734       fPIDSel->Fill(0);
 735         if (ptPart >= ptBinH[0] && ptPart < ptBinH[1])      fPIDSel_Bin1->Fill(0);
 736         else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fPIDSel_Bin2->Fill(0);
 737         else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fPIDSel_Bin3->Fill(0);
 738         else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fPIDSel_Bin4->Fill(0);
 739         else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fPIDSel_Bin5->Fill(0);
 740         }
 741
 742
 743    PostData(6,fOutputC);
 744
 745    //selection
 746    if((fSelected==1||fSelected==2||fSelected==3) && flagAccLim == 1) {
 747       // get daughter AOD tracks
 748       AliAODTrack *trk0 = (AliAODTrack*)dIn->GetDaughter(0);
 749       AliAODTrack *trk1 = (AliAODTrack*)dIn->GetDaughter(1);
 750       AliAODTrack *trk2 = (AliAODTrack*)dIn->GetDaughter(2);
 751       AliAODTrack *trk3 = (AliAODTrack*)dIn->GetDaughter(3);
 752       if(!trk0 || !trk1 || !trk2 || !trk3) {
 753         trk0=aodIn->GetTrack(trkIDtoEntry[dIn->GetProngID(0)]);
 754         trk1=aodIn->GetTrack(trkIDtoEntry[dIn->GetProngID(1)]);
 755         trk2=aodIn->GetTrack(trkIDtoEntry[dIn->GetProngID(2)]);
 756         trk3=aodIn->GetTrack(trkIDtoEntry[dIn->GetProngID(3)]);
 757       }
 758       printf("pt of positive track #1: %f\n",trk0->Pt());
 759       printf("pt of negative track #1: %f\n",trk1->Pt());
 760       printf("pt of positive track #2: %f\n",trk2->Pt());
 761       printf("pt of negative track #2: %f\n",trk3->Pt());
 762
 763       dIn->InvMassD0(fmassD0);
 764       dIn->InvMassD0bar(fmassD0bar);
 765
 766       //fill histos (combining selection form cuts & PID (with rho information))
 767       Int_t hypD01 = 0, hypD02 = 0, hypD0bar1 = 0, hypD0bar2 = 0;
 768       Int_t pid1 = 0, pid2 = 0, pidbar1 = 0, pidbar2 = 0;
 769       Int_t pidSelection = fCuts->IsSelectedFromPID(dIn, &pid1, &pid2, &pidbar1, &pidbar2);
 770
 771    //number of CANDIDATES (regardless of hypotheses) passing PPR + PID - PAY ATTENTION: hypoth. for PID and PPR may not be the same!
 772    if (pidSelection > 0) {
 773         fPIDSel->Fill(1);
 774         if (ptPart >= ptBinH[0] && ptPart < ptBinH[1])      fPIDSel_Bin1->Fill(1);
 775         else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fPIDSel_Bin2->Fill(1);
 776         else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fPIDSel_Bin3->Fill(1);
 777         else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fPIDSel_Bin4->Fill(1);
 778         else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fPIDSel_Bin5->Fill(1);
 779       }
 780
 781            //number of hypoteses accepted per candidate after PPR
 782            if(selD01+selD02+selD0bar1+selD0bar2 == 1) fMultipleHyps->Fill(0);
 783            if(selD01+selD02+selD0bar1+selD0bar2 == 2) fMultipleHyps->Fill(1);
 784            if(selD01+selD02+selD0bar1+selD0bar2 == 3) fMultipleHyps->Fill(2);
 785            if(selD01+selD02+selD0bar1+selD0bar2 == 4) fMultipleHyps->Fill(3);
 786
 787    //combine PPD + PID cuts
 788    if (selD01 == 1 && pid1==1) hypD01 = 1;
 789    if (selD02 == 1 && pid2==1) hypD02 = 1;
 790    if (selD0bar1 == 1 && pidbar1==1) hypD0bar1 = 1;
 791    if (selD0bar2 == 1 && pidbar2==1) hypD0bar2 = 1;
 792
 793    //number of CANDIDATES (regardless of hypotheses) passing PPR + PID - PAY ATTENTION: hypoth. for PID and PPR must match (at least one)!
 794    if (hypD01 == 1 || hypD02 == 1 || hypD0bar1 == 1 || hypD0bar2 == 1) {
 795         fPIDSel->Fill(2);
 796         if (ptPart >= ptBinH[0] && ptPart < ptBinH[1])      fPIDSel_Bin1->Fill(2);
 797         else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fPIDSel_Bin2->Fill(2);
 798         else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fPIDSel_Bin3->Fill(2);
 799         else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fPIDSel_Bin4->Fill(2);
 800         else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fPIDSel_Bin5->Fill(2);
 801       }
 802
 803            //number of hypoteses accepted per candidate after PPR and PID
 804            if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 1) fMultipleHyps->Fill(4);
 805            if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 2) fMultipleHyps->Fill(5);
 806            if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 3) fMultipleHyps->Fill(6);
 807            if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 4) fMultipleHyps->Fill(7);
 808
 809            //type of hypoteses accepted per candidate after PPR and PID
 810            if(hypD01+hypD02 == 1 && hypD0bar1+hypD0bar2 == 0) fMultipleHypsType->Fill(0);
 811            if(hypD01+hypD02 == 0 && hypD0bar1+hypD0bar2 == 1) fMultipleHypsType->Fill(1);
 812            if(hypD01+hypD02 == 2 && hypD0bar1+hypD0bar2 == 0) fMultipleHypsType->Fill(2);
 813            if(hypD01+hypD02 == 0 && hypD0bar1+hypD0bar2 == 2) fMultipleHypsType->Fill(3);
 814            if(hypD01+hypD02 == 1 && hypD0bar1+hypD0bar2 == 1) fMultipleHypsType->Fill(4);
 815            if(hypD01+hypD02 == 2 && hypD0bar1+hypD0bar2 == 1) fMultipleHypsType->Fill(5);
 816            if(hypD01+hypD02 == 1 && hypD0bar1+hypD0bar2 == 2) fMultipleHypsType->Fill(6);
 817            if(hypD01+hypD02 == 2 && hypD0bar1+hypD0bar2 == 2) fMultipleHypsType->Fill(7);
 818
 819
 820 //
 821 // This section is auxiliary: (multiple hypotheses histos)
 822 //
 823  if (hypD01+hypD02+hypD0bar1+hypD0bar2 == 2 || hypD01+hypD02+hypD0bar1+hypD0bar2 == 3 || hypD01+hypD02+hypD0bar1+hypD0bar2 == 4)
 824  {
 825
 826    if (ptPart > ptBinH[0]) {
 827       // D01 hyp.
 828       if(hypD01==1) {
 829         if (ptPart < ptBinH[1])                             fhInvMassMultipleOnly_Bin1->Fill(fmassD0[0]);
 830         else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fhInvMassMultipleOnly_Bin2->Fill(fmassD0[0]);
 831         else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fhInvMassMultipleOnly_Bin3->Fill(fmassD0[0]);
 832         else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fhInvMassMultipleOnly_Bin4->Fill(fmassD0[0]);
 833         else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fhInvMassMultipleOnly_Bin5->Fill(fmassD0[0]);
 834        }
 835       // D02 hyp.
 836       if(hypD02==1) {
 837         if (ptPart < ptBinH[1])                             fhInvMassMultipleOnly_Bin1->Fill(fmassD0[1]);
 838         else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fhInvMassMultipleOnly_Bin2->Fill(fmassD0[1]);
 839         else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fhInvMassMultipleOnly_Bin3->Fill(fmassD0[1]);
 840         else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fhInvMassMultipleOnly_Bin4->Fill(fmassD0[1]);
 841         else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fhInvMassMultipleOnly_Bin5->Fill(fmassD0[1]);
 842        }
 843       // D0bar1 hyp.
 844       if(hypD0bar1==1) {
 845         if (ptPart < ptBinH[1])                             fhInvMassMultipleOnly_Bin1->Fill(fmassD0bar[0]);
 846         else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fhInvMassMultipleOnly_Bin2->Fill(fmassD0bar[0]);
 847         else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fhInvMassMultipleOnly_Bin3->Fill(fmassD0bar[0]);
 848         else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fhInvMassMultipleOnly_Bin4->Fill(fmassD0bar[0]);
 849         else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fhInvMassMultipleOnly_Bin5->Fill(fmassD0bar[0]);
 850        }
 851       // D0bar2 hyp.
 852       if(hypD0bar2==1) {
 853         if (ptPart < ptBinH[1])                             fhInvMassMultipleOnly_Bin1->Fill(fmassD0bar[1]);
 854         else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) fhInvMassMultipleOnly_Bin2->Fill(fmassD0bar[1]);
 855         else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) fhInvMassMultipleOnly_Bin3->Fill(fmassD0bar[1]);
 856         else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) fhInvMassMultipleOnly_Bin4->Fill(fmassD0bar[1]);
 857         else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5]) fhInvMassMultipleOnly_Bin5->Fill(fmassD0bar[1]);
 858        }
 859      }
 860
 861    }
 862 //
 863 //end of auxiliary section
 864 //
 865
 866
 867    //All histos are filled if Pt of candidate is greater than minimum of first bin (in this way: bin1+bin2+...binN = whole)
 868    if (ptPart > ptBinH[0]) {
 869
 870       // D01 hyp.
 871       if(hypD01==1) {
 872         fPtSel->Fill(ptPart);
 873         fScatterP4PID->Fill(trk1->Pt(),trk3->Pt());
 874         if (ptPart < ptBinH[1])                          {fhInvMassD0Sum_10Mev_Bin1->Fill(fmassD0[0]);
 875                                                           fhInvMassD0Sum_5Mev_Bin1->Fill(fmassD0[0]);
 876                                                           fhInvMassSumAll_10Mev_Bin1->Fill(fmassD0[0]);
 877                                                           fhInvMassSumAll_5Mev_Bin1->Fill(fmassD0[0]);}
 878         else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) {fhInvMassD0Sum_10Mev_Bin2->Fill(fmassD0[0]);
 879                                                              fhInvMassD0Sum_5Mev_Bin2->Fill(fmassD0[0]);
 880                                                              fhInvMassSumAll_10Mev_Bin2->Fill(fmassD0[0]);
 881                                                              fhInvMassSumAll_5Mev_Bin2->Fill(fmassD0[0]);}
 882         else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) {fhInvMassD0Sum_10Mev_Bin3->Fill(fmassD0[0]);
 883                                                              fhInvMassD0Sum_5Mev_Bin3->Fill(fmassD0[0]);
 884                                                              fhInvMassSumAll_10Mev_Bin3->Fill(fmassD0[0]);
 885                                                              fhInvMassSumAll_5Mev_Bin3->Fill(fmassD0[0]);}
 886         else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) {fhInvMassD0Sum_10Mev_Bin4->Fill(fmassD0[0]);
 887                                                              fhInvMassD0Sum_5Mev_Bin4->Fill(fmassD0[0]);
 888                                                              fhInvMassSumAll_10Mev_Bin4->Fill(fmassD0[0]);
 889                                                              fhInvMassSumAll_5Mev_Bin4->Fill(fmassD0[0]);}
 890         else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5])     {fhInvMassD0Sum_10Mev_Bin5->Fill(fmassD0[0]);
 891                                                                 fhInvMassD0Sum_5Mev_Bin5->Fill(fmassD0[0]);
 892                                                                 fhInvMassSumAll_10Mev_Bin5->Fill(fmassD0[0]);
 893                                                                 fhInvMassSumAll_5Mev_Bin5->Fill(fmassD0[0]);}
 894      }
 895       // D02 hyp.
 896       if(hypD02==1) {
 897         fPtSel->Fill(ptPart);
 898         fScatterP4PID->Fill(trk3->Pt(),trk1->Pt());
 899         if (ptPart < ptBinH[1])                          {fhInvMassD0Sum_10Mev_Bin1->Fill(fmassD0[1]);
 900                                                           fhInvMassD0Sum_5Mev_Bin1->Fill(fmassD0[1]);
 901                                                           fhInvMassSumAll_10Mev_Bin1->Fill(fmassD0[1]);
 902                                                           fhInvMassSumAll_5Mev_Bin1->Fill(fmassD0[1]);}
 903         else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) {fhInvMassD0Sum_10Mev_Bin2->Fill(fmassD0[1]);
 904                                                              fhInvMassD0Sum_5Mev_Bin2->Fill(fmassD0[1]);
 905                                                              fhInvMassSumAll_10Mev_Bin2->Fill(fmassD0[1]);
 906                                                              fhInvMassSumAll_5Mev_Bin2->Fill(fmassD0[1]);}
 907         else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) {fhInvMassD0Sum_10Mev_Bin3->Fill(fmassD0[1]);
 908                                                              fhInvMassD0Sum_5Mev_Bin3->Fill(fmassD0[1]);
 909                                                              fhInvMassSumAll_10Mev_Bin3->Fill(fmassD0[1]);
 910                                                              fhInvMassSumAll_5Mev_Bin3->Fill(fmassD0[1]);}
 911         else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) {fhInvMassD0Sum_10Mev_Bin4->Fill(fmassD0[1]);
 912                                                              fhInvMassD0Sum_5Mev_Bin4->Fill(fmassD0[1]);
 913                                                              fhInvMassSumAll_10Mev_Bin4->Fill(fmassD0[1]);
 914                                                              fhInvMassSumAll_5Mev_Bin4->Fill(fmassD0[1]);}
 915         else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5])     {fhInvMassD0Sum_10Mev_Bin5->Fill(fmassD0[1]);
 916                                                                  fhInvMassD0Sum_5Mev_Bin5->Fill(fmassD0[1]);
 917                                                                  fhInvMassSumAll_10Mev_Bin5->Fill(fmassD0[1]);
 918                                                                  fhInvMassSumAll_5Mev_Bin5->Fill(fmassD0[1]);}
 919      }
 920       // D0bar1 hyp.
 921       if(hypD0bar1==1) {
 922         fPtSel->Fill(ptPart);
 923         fScatterP4PID->Fill(trk0->Pt(),trk2->Pt());
 924         if (ptPart < ptBinH[1])                          {fhInvMassD0barSum_10Mev_Bin1->Fill(fmassD0bar[0]);
 925                                                           fhInvMassD0barSum_5Mev_Bin1->Fill(fmassD0bar[0]);
 926                                                           fhInvMassSumAll_10Mev_Bin1->Fill(fmassD0bar[0]);
 927                                                           fhInvMassSumAll_5Mev_Bin1->Fill(fmassD0bar[0]);}
 928         else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) {fhInvMassD0barSum_10Mev_Bin2->Fill(fmassD0bar[0]);
 929                                                              fhInvMassD0barSum_5Mev_Bin2->Fill(fmassD0bar[0]);
 930                                                              fhInvMassSumAll_10Mev_Bin2->Fill(fmassD0bar[0]);
 931                                                              fhInvMassSumAll_5Mev_Bin2->Fill(fmassD0bar[0]);}
 932         else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) {fhInvMassD0barSum_10Mev_Bin3->Fill(fmassD0bar[0]);
 933                                                              fhInvMassD0barSum_5Mev_Bin3->Fill(fmassD0bar[0]);
 934                                                              fhInvMassSumAll_10Mev_Bin3->Fill(fmassD0bar[0]);
 935                                                              fhInvMassSumAll_5Mev_Bin3->Fill(fmassD0bar[0]);}
 936         else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) {fhInvMassD0barSum_10Mev_Bin4->Fill(fmassD0bar[0]);
 937                                                              fhInvMassD0barSum_5Mev_Bin4->Fill(fmassD0bar[0]);
 938                                                              fhInvMassSumAll_10Mev_Bin4->Fill(fmassD0bar[0]);
 939                                                              fhInvMassSumAll_5Mev_Bin4->Fill(fmassD0bar[0]);}
 940         else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5])     {fhInvMassD0barSum_10Mev_Bin5->Fill(fmassD0bar[0]);
 941                                                                 fhInvMassD0barSum_5Mev_Bin5->Fill(fmassD0bar[0]);
 942                                                                 fhInvMassSumAll_10Mev_Bin5->Fill(fmassD0bar[0]);
 943                                                                 fhInvMassSumAll_5Mev_Bin5->Fill(fmassD0bar[0]);}
 944      }
 945       // D0bar2 hyp.
 946       if(hypD0bar2==1) {
 947         fPtSel->Fill(ptPart);
 948         fScatterP4PID->Fill(trk2->Pt(),trk0->Pt());
 949         if (ptPart < ptBinH[1])                          {fhInvMassD0barSum_10Mev_Bin1->Fill(fmassD0bar[1]);
 950                                                           fhInvMassD0barSum_5Mev_Bin1->Fill(fmassD0bar[1]);
 951                                                           fhInvMassSumAll_10Mev_Bin1->Fill(fmassD0bar[1]);
 952                                                           fhInvMassSumAll_5Mev_Bin1->Fill(fmassD0bar[1]);}
 953         else if (ptPart >= ptBinH[1] && ptPart < ptBinH[2]) {fhInvMassD0barSum_10Mev_Bin2->Fill(fmassD0bar[1]);
 954                                                              fhInvMassD0barSum_5Mev_Bin2->Fill(fmassD0bar[1]);
 955                                                              fhInvMassSumAll_10Mev_Bin2->Fill(fmassD0bar[1]);
 956                                                              fhInvMassSumAll_5Mev_Bin2->Fill(fmassD0bar[1]);}
 957         else if (ptPart >= ptBinH[2] && ptPart < ptBinH[3]) {fhInvMassD0barSum_10Mev_Bin3->Fill(fmassD0bar[1]);
 958                                                              fhInvMassD0barSum_5Mev_Bin3->Fill(fmassD0bar[1]);
 959                                                              fhInvMassSumAll_10Mev_Bin3->Fill(fmassD0bar[1]);
 960                                                              fhInvMassSumAll_5Mev_Bin3->Fill(fmassD0bar[1]);}
 961         else if (ptPart >= ptBinH[3] && ptPart < ptBinH[4]) {fhInvMassD0barSum_10Mev_Bin4->Fill(fmassD0bar[1]);
 962                                                              fhInvMassD0barSum_5Mev_Bin4->Fill(fmassD0bar[1]);
 963                                                              fhInvMassSumAll_10Mev_Bin4->Fill(fmassD0bar[1]);
 964                                                              fhInvMassSumAll_5Mev_Bin4->Fill(fmassD0bar[1]);}
 965         else if (ptPart >= ptBinH[4] && ptPart < ptBinH[5])     {fhInvMassD0barSum_10Mev_Bin5->Fill(fmassD0bar[1]);
 966                                                                 fhInvMassD0barSum_5Mev_Bin5->Fill(fmassD0bar[1]);
 967                                                                 fhInvMassSumAll_10Mev_Bin5->Fill(fmassD0bar[1]);
 968                                                                 fhInvMassSumAll_5Mev_Bin5->Fill(fmassD0bar[1]);}
 969      }
 970    }
 971
 972       PostData(1,fOutput);
 973       PostData(2,fOutput2);
 974       PostData(3,fOutput3);
 975       PostData(4,fOutput4);
 976       PostData(5,fOutput5);
 977
 978       // HERE ONE COULD RECALCULATE THE VERTEX USING THE KF PACKAGE
 979
 980       // clone candidate for output AOD
 981       if(hypD01||hypD02||hypD0bar1||hypD0bar2) {
 982       AliAODVertex *v = new(verticesHFRef[iOutVerticesHF++])
 983         AliAODVertex(*(dIn->GetSecondaryVtx()));
 984       AliAODRecoDecayHF4Prong *dOut=new(aodCharm4ProngRef[iOutCharm4Prong++])
 985         AliAODRecoDecayHF4Prong(*dIn);
 986       dOut->SetSecondaryVtx(v);
 987       dOut->SetOwnPrimaryVtx((AliAODVertex*)((dIn->GetOwnPrimaryVtx())->Clone()));
 988       v->SetParent(dOut);
 989       }
 990     } //end of selection loop (starts with fSelected == 1, 2 or 3)
 991     if(unsetvtx) dIn->UnsetOwnPrimaryVtx();
 992   } // end loop on D0->K3pi
 993
 994   printf("Number of selected D0->K3pi: %d\n",iOutCharm4Prong);
 995
 996   return;
 997 }
 998
 999 //________________________________________________________________________
1000 void AliAnalysisTaskSESelectHF4Prong::Terminate(Option_t */*option*/)
1001 {
1002   // Terminate analysis
1003   //
1004   if(fDebug > 1) printf("AnalysisTaskSESelectHF4Prong: Terminate() \n");
1005
1006   fOutput = dynamic_cast<TList*> (GetOutputData(1));
1007   if (!fOutput) {
1008     printf("ERROR: fOutput not available\n");
1009     return;
1010   }
1011   fOutput2 = dynamic_cast<TList*> (GetOutputData(2));
1012   if (!fOutput2) {
1013     printf("ERROR: fOutput not available\n");
1014     return;
1015   }
1016   fOutput3 = dynamic_cast<TList*> (GetOutputData(3));
1017   if (!fOutput3) {
1018     printf("ERROR: fOutput not available\n");
1019     return;
1020   }
1021   fOutput4 = dynamic_cast<TList*> (GetOutputData(4));
1022   if (!fOutput4) {
1023     printf("ERROR: fOutput not available\n");
1024     return;
1025   }
1026   if (!fOutput5) {
1027     printf("ERROR: fOutput not available\n");
1028     return;
1029   }
1030   fOutputC = dynamic_cast<TList*> (GetOutputData(5));
1031   if (!fOutputC) {
1032     printf("ERROR: fOutputC not available\n");
1033     return;
1034   }
1035
1036   fhInvMassD0Sum_10Mev_Bin1 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_10Mev_Bin1"));
1037   fhInvMassD0barSum_10Mev_Bin1 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_10Mev_Bin1"));
1038   fhInvMassSumAll_10Mev_Bin1 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_10Mev_Bin1"));
1039   fhInvMassD0Sum_5Mev_Bin1 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_5Mev_Bin1"));
1040   fhInvMassD0barSum_5Mev_Bin1 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_5Mev_Bin1"));
1041   fhInvMassSumAll_5Mev_Bin1 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_5Mev_Bin1"));
1042   fhInvMassD0Sum_10Mev_Bin2 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_10Mev_Bin2"));
1043   fhInvMassD0barSum_10Mev_Bin2 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_10Mev_Bin2"));
1044   fhInvMassSumAll_10Mev_Bin2 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_10Mev_Bin2"));
1045   fhInvMassD0Sum_5Mev_Bin2 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_5Mev_Bin2"));
1046   fhInvMassD0barSum_5Mev_Bin2 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_5Mev_Bin2"));
1047   fhInvMassSumAll_5Mev_Bin2 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_5Mev_Bin2"));
1048   fhInvMassD0Sum_10Mev_Bin3 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_10Mev_Bin3"));
1049   fhInvMassD0barSum_10Mev_Bin3 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_10Mev_Bin3"));
1050   fhInvMassSumAll_10Mev_Bin3 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_10Mev_Bin3"));
1051   fhInvMassD0Sum_5Mev_Bin3 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_5Mev_Bin3"));
1052   fhInvMassD0barSum_5Mev_Bin3 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_5Mev_Bin3"));
1053   fhInvMassSumAll_5Mev_Bin3 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_5Mev_Bin3"));
1054   fhInvMassD0Sum_10Mev_Bin4 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_10Mev_Bin4"));
1055   fhInvMassD0barSum_10Mev_Bin4 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_10Mev_Bin4"));
1056   fhInvMassSumAll_10Mev_Bin4 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_10Mev_Bin4"));
1057   fhInvMassD0Sum_5Mev_Bin4 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_5Mev_Bin4"));
1058   fhInvMassD0barSum_5Mev_Bin4 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_5Mev_Bin4"));
1059   fhInvMassSumAll_5Mev_Bin4 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_5Mev_Bin4"));
1060   fhInvMassD0Sum_10Mev_Bin5 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_10Mev_Bin5"));
1061   fhInvMassD0barSum_10Mev_Bin5 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_10Mev_Bin5"));
1062   fhInvMassSumAll_10Mev_Bin5 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_10Mev_Bin5"));
1063   fhInvMassD0Sum_5Mev_Bin5 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum_5Mev_Bin5"));
1064   fhInvMassD0barSum_5Mev_Bin5 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum_5Mev_Bin5"));
1065   fhInvMassSumAll_5Mev_Bin5 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll_5Mev_Bin5"));
1066
1067   fScatterP4PID = dynamic_cast<TH2F*>(fOutput->FindObject("fScatterP4PID"));
1068   fPtVsY = dynamic_cast<TH2F*>(fOutputC->FindObject("fPtVsY"));
1069   fPtVsYAll = dynamic_cast<TH2F*>(fOutputC->FindObject("fPtVsYAll"));
1070
1071   fEventCounter = dynamic_cast<TH1F*>(fOutputC->FindObject("fEventCounter"));
1072
1073   fCutDCA = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutDCA"));
1074   fCutDCA3 = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutDCA3"));
1075   fCutDCA2 = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutDCA2"));
1076   fCutDCA5 = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutDCA5"));
1077   fCutVertexDist2 = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutVertexDist2"));
1078   fCutVertexDist3 = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutVertexDist3"));
1079   fCutVertexDist4 = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutVertexDist4"));
1080   fCutCosinePoint = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutCosinePoint"));
1081
1082   fCutPt = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutPt"));
1083   fCutY = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutY"));
1084   fPIDSel = dynamic_cast<TH1F*>(fOutputC->FindObject("fPIDSel"));
1085   fPIDSel_Bin1 = dynamic_cast<TH1F*>(fOutputC->FindObject("fPIDSel_Bin1"));
1086   fPIDSel_Bin2 = dynamic_cast<TH1F*>(fOutputC->FindObject("fPIDSel_Bin2"));
1087   fPIDSel_Bin3 = dynamic_cast<TH1F*>(fOutputC->FindObject("fPIDSel_Bin3"));
1088   fPIDSel_Bin4 = dynamic_cast<TH1F*>(fOutputC->FindObject("fPIDSel_Bin4"));
1089   fPIDSel_Bin5 = dynamic_cast<TH1F*>(fOutputC->FindObject("fPIDSel_Bin5"));
1090 }
1091