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 fSelected(0),
  62 fMCTruth(0),
  63 fOutput(0),
  64 fOutput2(0),
  65 fOutput3(0),
  66 fOutput4(0),
  67 fOutput5(0),
  68 fOutputC(0),
  69 fhInvMassD0Sum10MevBin1(0),
  70 fhInvMassD0barSum10MevBin1(0),
  71 fhInvMassSumAll10MevBin1(0),
  72 fhInvMassD0Sum5MevBin1(0),
  73 fhInvMassD0barSum5MevBin1(0),
  74 fhInvMassSumAll5MevBin1(0),
  75 fhInvMassD0Sum10MevBin2(0),
  76 fhInvMassD0barSum10MevBin2(0),
  77 fhInvMassSumAll10MevBin2(0),
  78 fhInvMassD0Sum5MevBin2(0),
  79 fhInvMassD0barSum5MevBin2(0),
  80 fhInvMassSumAll5MevBin2(0),
  81 fhInvMassD0Sum10MevBin3(0),
  82 fhInvMassD0barSum10MevBin3(0),
  83 fhInvMassSumAll10MevBin3(0),
  84 fhInvMassD0Sum5MevBin3(0),
  85 fhInvMassD0barSum5MevBin3(0),
  86 fhInvMassSumAll5MevBin3(0),
  87 fhInvMassD0Sum10MevBin4(0),
  88 fhInvMassD0barSum10MevBin4(0),
  89 fhInvMassSumAll10MevBin4(0),
  90 fhInvMassD0Sum5MevBin4(0),
  91 fhInvMassD0barSum5MevBin4(0),
  92 fhInvMassSumAll5MevBin4(0),
  93 fhInvMassD0Sum10MevBin5(0),
  94 fhInvMassD0barSum10MevBin5(0),
  95 fhInvMassSumAll10MevBin5(0),
  96 fhInvMassD0Sum5MevBin5(0),
  97 fhInvMassD0barSum5MevBin5(0),
  98 fhInvMassSumAll5MevBin5(0),
  99 fhReflBin1(0),
 100 fhReflBin2(0),
 101 fhReflBin3(0),
 102 fhReflBin4(0),
 103 fhReflBin5(0),
 104 fhReflD0Bin1(0),
 105 fhReflD0Bin2(0),
 106 fhReflD0Bin3(0),
 107 fhReflD0Bin4(0),
 108 fhReflD0Bin5(0),
 109 fhReflD0barBin1(0),
 110 fhReflD0barBin2(0),
 111 fhReflD0barBin3(0),
 112 fhReflD0barBin4(0),
 113 fhReflD0barBin5(0),
 114 fScatterP4PID(0),
 115 fPtVsY(0),
 116 fPtVsYAll(0),
 117 fEventCounter(0),
 118 fCutDCA(0),
 119 fCutDCA3(0),
 120 fCutDCA2(0),
 121 fCutDCA5(0),
 122 fCutVertexDist2(0),
 123 fCutVertexDist3(0),
 124 fCutVertexDist4(0),
 125 fCutCosinePoint(0),
 126 fCutPt(0),
 127 fCutY(0),
 128 fPIDSel(0),
 129 fPIDSelBin1(0),
 130 fPIDSelBin2(0),
 131 fPIDSelBin3(0),
 132 fPIDSelBin4(0),
 133 fPIDSelBin5(0),
 134 fMultipleHyps(0),
 135 fMultipleHypsType(0),
 136 fPtSel(0),
 137 fCuts(0)
 138 {
 139   // Default constructor
 140 }
 141
 142 //_______________________________________________________________________
 143 AliAnalysisTaskSESelectHF4Prong::AliAnalysisTaskSESelectHF4Prong(const char *name,AliRDHFCutsD0toKpipipi* cuts):
 144 AliAnalysisTaskSE(name),
 145 fVerticesHFTClArr(0),
 146 fCharm4ProngTClArr(0),
 147 fSelected(0),
 148 fMCTruth(0),
 149 fOutput(0),
 150 fOutput2(0),
 151 fOutput3(0),
 152 fOutput4(0),
 153 fOutput5(0),
 154 fOutputC(0),
 155 fhInvMassD0Sum10MevBin1(0),
 156 fhInvMassD0barSum10MevBin1(0),
 157 fhInvMassSumAll10MevBin1(0),
 158 fhInvMassD0Sum5MevBin1(0),
 159 fhInvMassD0barSum5MevBin1(0),
 160 fhInvMassSumAll5MevBin1(0),
 161 fhInvMassD0Sum10MevBin2(0),
 162 fhInvMassD0barSum10MevBin2(0),
 163 fhInvMassSumAll10MevBin2(0),
 164 fhInvMassD0Sum5MevBin2(0),
 165 fhInvMassD0barSum5MevBin2(0),
 166 fhInvMassSumAll5MevBin2(0),
 167 fhInvMassD0Sum10MevBin3(0),
 168 fhInvMassD0barSum10MevBin3(0),
 169 fhInvMassSumAll10MevBin3(0),
 170 fhInvMassD0Sum5MevBin3(0),
 171 fhInvMassD0barSum5MevBin3(0),
 172 fhInvMassSumAll5MevBin3(0),
 173 fhInvMassD0Sum10MevBin4(0),
 174 fhInvMassD0barSum10MevBin4(0),
 175 fhInvMassSumAll10MevBin4(0),
 176 fhInvMassD0Sum5MevBin4(0),
 177 fhInvMassD0barSum5MevBin4(0),
 178 fhInvMassSumAll5MevBin4(0),
 179 fhInvMassD0Sum10MevBin5(0),
 180 fhInvMassD0barSum10MevBin5(0),
 181 fhInvMassSumAll10MevBin5(0),
 182 fhInvMassD0Sum5MevBin5(0),
 183 fhInvMassD0barSum5MevBin5(0),
 184 fhInvMassSumAll5MevBin5(0),
 185 fhReflBin1(0),
 186 fhReflBin2(0),
 187 fhReflBin3(0),
 188 fhReflBin4(0),
 189 fhReflBin5(0),
 190 fhReflD0Bin1(0),
 191 fhReflD0Bin2(0),
 192 fhReflD0Bin3(0),
 193 fhReflD0Bin4(0),
 194 fhReflD0Bin5(0),
 195 fhReflD0barBin1(0),
 196 fhReflD0barBin2(0),
 197 fhReflD0barBin3(0),
 198 fhReflD0barBin4(0),
 199 fhReflD0barBin5(0),
 200 fScatterP4PID(0),
 201 fPtVsY(0),
 202 fPtVsYAll(0),
 203 fEventCounter(0),
 204 fCutDCA(0),
 205 fCutDCA3(0),
 206 fCutDCA2(0),
 207 fCutDCA5(0),
 208 fCutVertexDist2(0),
 209 fCutVertexDist3(0),
 210 fCutVertexDist4(0),
 211 fCutCosinePoint(0),
 212 fCutPt(0),
 213 fCutY(0),
 214 fPIDSel(0),
 215 fPIDSelBin1(0),
 216 fPIDSelBin2(0),
 217 fPIDSelBin3(0),
 218 fPIDSelBin4(0),
 219 fPIDSelBin5(0),
 220 fMultipleHyps(0),
 221 fMultipleHypsType(0),
 222 fPtSel(0),
 223 fCuts(0)
 224 {
 225   // Standard constructor
 226
 227   fCuts=cuts;
 228
 229   // Input slot #0 works with an Ntuple
 230 //  DefineInput(0, TTree::Class());
 231
 232   // Output slot #0 writes into a TTree container
 233   // Output slots #1-6 writes into a TList container
 234   DefineOutput(0, TTree::Class()); //default
 235   DefineOutput(1, TList::Class()); //histos inv. mass bin1
 236   DefineOutput(2, TList::Class()); //histos inv. mass bin2
 237   DefineOutput(3, TList::Class()); //histos inv. mass bin3
 238   DefineOutput(4, TList::Class()); //histos inv. mass bin4
 239   DefineOutput(5, TList::Class()); //histos inv. mass bin5
 240   DefineOutput(6, TList::Class()); //histos of cuts
 241   DefineOutput(7, AliRDHFCutsD0toKpipipi::Class()); //cuts
 242 }
 243
 244 //________________________________________________________________________
 245 AliAnalysisTaskSESelectHF4Prong::~AliAnalysisTaskSESelectHF4Prong()
 246 {
 247   // Destructor
 248
 249   if (fOutput) {
 250     delete fOutput;
 251     fOutput = 0;
 252   }
 253   if (fOutput2) {
 254     delete fOutput2;
 255     fOutput2 = 0;
 256   }
 257   if (fOutput3) {
 258     delete fOutput3;
 259     fOutput3 = 0;
 260   }
 261   if (fOutput4) {
 262     delete fOutput4;
 263     fOutput4 = 0;
 264   }
 265   if (fOutput5) {
 266     delete fOutput5;
 267     fOutput5 = 0;
 268   }
 269   if (fOutputC) {
 270     delete fOutputC;
 271     fOutputC = 0;
 272   }
 273   if (fCuts) {
 274     delete fCuts;
 275     fCuts = 0;
 276   }
 277
 278 }
 279
 280 //________________________________________________________________________
 281 void AliAnalysisTaskSESelectHF4Prong::Init()
 282 {
 283   // Initialization
 284
 285   if(fDebug > 1) printf("AnalysisTaskSESelectHF4Prong::Init() \n");
 286   PrintPtBinHandMCFlag();
 287
 288   return;
 289 }
 290
 291 //________________________________________________________________________
 292 void AliAnalysisTaskSESelectHF4Prong::UserCreateOutputObjects()
 293 {
 294   // Create the output container
 295   //
 296   if(fDebug > 1) printf("AnalysisTaskSESelectHF4Prong::UserCreateOutputObjects() \n");
 297
 298   if(fDebug > 1) PrintPtBinHandMCFlag();
 299
 300   fVerticesHFTClArr = new TClonesArray("AliAODVertex", 0);
 301   fVerticesHFTClArr->SetName("VerticesHF");
 302   AddAODBranch("TClonesArray", &fVerticesHFTClArr);
 303
 304   fCharm4ProngTClArr = new TClonesArray("AliAODRecoDecayHF4Prong", 0);
 305   fCharm4ProngTClArr->SetName("Charm4Prong");
 306   AddAODBranch("TClonesArray", &fCharm4ProngTClArr);
 307
 308   fOutput = new TList();
 309   fOutput->SetOwner();
 310
 311   fOutput2 = new TList();
 312   fOutput2->SetOwner();
 313
 314   fOutput3 = new TList();
 315   fOutput3->SetOwner();
 316
 317   fOutput4 = new TList();
 318   fOutput4->SetOwner();
 319
 320   fOutput5 = new TList();
 321   fOutput5->SetOwner();
 322
 323   fOutputC = new TList();
 324   fOutputC->SetOwner();
 325
 326   fhInvMassD0Sum10MevBin1 = new TH1F("fhInvMassD0Sum10MevBin1", "D0 invariant mass Bin1 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
 327   fhInvMassD0Sum10MevBin1->Sumw2(); //Create structure to store sum of squares of weights
 328   fhInvMassD0Sum10MevBin1->SetMinimum(0);
 329   fOutput->Add(fhInvMassD0Sum10MevBin1);
 330
 331   fhInvMassD0barSum10MevBin1 = new TH1F("fhInvMassD0barSum10MevBin1", "D0bar invariant mass Bin1 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
 332   fhInvMassD0barSum10MevBin1->Sumw2(); //Create structure to store sum of squares of weights
 333   fhInvMassD0barSum10MevBin1->SetMinimum(0);
 334   fOutput->Add(fhInvMassD0barSum10MevBin1);
 335
 336   fhInvMassSumAll10MevBin1 = new TH1F("fhInvMassSumAll10MevBin1", "D0/D0bar invariant mass Bin1 (good hyps); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
 337   fhInvMassSumAll10MevBin1->Sumw2(); //Create structure to store sum of squares of weights
 338   fhInvMassSumAll10MevBin1->SetMinimum(0);
 339   fOutput->Add(fhInvMassSumAll10MevBin1);
 340
 341   fhInvMassD0Sum5MevBin1 = new TH1F("fhInvMassD0Sum5MevBin1", "D0 invariant mass Bin1 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
 342   fhInvMassD0Sum5MevBin1->Sumw2(); //Create structure to store sum of squares of weights
 343   fhInvMassD0Sum5MevBin1->SetMinimum(0);
 344   fOutput->Add(fhInvMassD0Sum5MevBin1);
 345
 346   fhInvMassD0barSum5MevBin1 = new TH1F("fhInvMassD0barSum5MevBin1", "D0bar invariant mass Bin1 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
 347   fhInvMassD0barSum5MevBin1->Sumw2(); //Create structure to store sum of squares of weights
 348   fhInvMassD0barSum5MevBin1->SetMinimum(0);
 349   fOutput->Add(fhInvMassD0barSum5MevBin1);
 350
 351   fhInvMassSumAll5MevBin1 = new TH1F("fhInvMassSumAll5MevBin1", "D0/D0bar invariant mass Bin1 (good hyps); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
 352   fhInvMassSumAll5MevBin1->Sumw2(); //Create structure to store sum of squares of weights
 353   fhInvMassSumAll5MevBin1->SetMinimum(0);
 354   fOutput->Add(fhInvMassSumAll5MevBin1);
 355
 356   fhInvMassD0Sum10MevBin2 = new TH1F("fhInvMassD0Sum10MevBin2", "D0 invariant mass Bin2 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
 357   fhInvMassD0Sum10MevBin2->Sumw2(); //Create structure to store sum of squares of weights
 358   fhInvMassD0Sum10MevBin2->SetMinimum(0);
 359   fOutput2->Add(fhInvMassD0Sum10MevBin2);
 360
 361   fhInvMassD0barSum10MevBin2 = new TH1F("fhInvMassD0barSum10MevBin2", "D0bar invariant mass Bin2 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
 362   fhInvMassD0barSum10MevBin2->Sumw2(); //Create structure to store sum of squares of weights
 363   fhInvMassD0barSum10MevBin2->SetMinimum(0);
 364   fOutput2->Add(fhInvMassD0barSum10MevBin2);
 365
 366   fhInvMassSumAll10MevBin2 = new TH1F("fhInvMassSumAll10MevBin2", "D0/D0bar invariant mass Bin2 (good hyps); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
 367   fhInvMassSumAll10MevBin2->Sumw2(); //Create structure to store sum of squares of weights
 368   fhInvMassSumAll10MevBin2->SetMinimum(0);
 369   fOutput2->Add(fhInvMassSumAll10MevBin2);
 370
 371   fhInvMassD0Sum5MevBin2 = new TH1F("fhInvMassD0Sum5MevBin2", "D0 invariant mass Bin2 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
 372   fhInvMassD0Sum5MevBin2->Sumw2(); //Create structure to store sum of squares of weights
 373   fhInvMassD0Sum5MevBin2->SetMinimum(0);
 374   fOutput2->Add(fhInvMassD0Sum5MevBin2);
 375
 376   fhInvMassD0barSum5MevBin2 = new TH1F("fhInvMassD0barSum5MevBin2", "D0bar invariant mass Bin2 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
 377   fhInvMassD0barSum5MevBin2->Sumw2(); //Create structure to store sum of squares of weights
 378   fhInvMassD0barSum5MevBin2->SetMinimum(0);
 379   fOutput2->Add(fhInvMassD0barSum5MevBin2);
 380
 381   fhInvMassSumAll5MevBin2 = new TH1F("fhInvMassSumAll5MevBin2", "D0/D0bar invariant mass Bin2 (good hyps); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
 382   fhInvMassSumAll5MevBin2->Sumw2(); //Create structure to store sum of squares of weights
 383   fhInvMassSumAll5MevBin2->SetMinimum(0);
 384   fOutput2->Add(fhInvMassSumAll5MevBin2);
 385
 386   fhInvMassD0Sum10MevBin3 = new TH1F("fhInvMassD0Sum10MevBin3", "D0 invariant mass Bin3 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
 387   fhInvMassD0Sum10MevBin3->Sumw2(); //Create structure to store sum of squares of weights
 388   fhInvMassD0Sum10MevBin3->SetMinimum(0);
 389   fOutput3->Add(fhInvMassD0Sum10MevBin3);
 390
 391   fhInvMassD0barSum10MevBin3 = new TH1F("fhInvMassD0barSum10MevBin3", "D0bar invariant mass Bin3 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
 392   fhInvMassD0barSum10MevBin3->Sumw2(); //Create structure to store sum of squares of weights
 393   fhInvMassD0barSum10MevBin3->SetMinimum(0);
 394   fOutput3->Add(fhInvMassD0barSum10MevBin3);
 395
 396   fhInvMassSumAll10MevBin3 = new TH1F("fhInvMassSumAll10MevBin3", "D0/D0bar invariant mass Bin3 (good hyps); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
 397   fhInvMassSumAll10MevBin3->Sumw2(); //Create structure to store sum of squares of weights
 398   fhInvMassSumAll10MevBin3->SetMinimum(0);
 399   fOutput3->Add(fhInvMassSumAll10MevBin3);
 400
 401   fhInvMassD0Sum5MevBin3 = new TH1F("fhInvMassD0Sum5MevBin3", "D0 invariant mass Bin3 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
 402   fhInvMassD0Sum5MevBin3->Sumw2(); //Create structure to store sum of squares of weights
 403   fhInvMassD0Sum5MevBin3->SetMinimum(0);
 404   fOutput3->Add(fhInvMassD0Sum5MevBin3);
 405
 406   fhInvMassD0barSum5MevBin3 = new TH1F("fhInvMassD0barSum5MevBin3", "D0bar invariant mass Bin3 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
 407   fhInvMassD0barSum5MevBin3->Sumw2(); //Create structure to store sum of squares of weights
 408   fhInvMassD0barSum5MevBin3->SetMinimum(0);
 409   fOutput3->Add(fhInvMassD0barSum5MevBin3);
 410
 411   fhInvMassSumAll5MevBin3 = new TH1F("fhInvMassSumAll5MevBin3", "D0/D0bar invariant mass Bin3 (good hyps); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
 412   fhInvMassSumAll5MevBin3->Sumw2(); //Create structure to store sum of squares of weights
 413   fhInvMassSumAll5MevBin3->SetMinimum(0);
 414   fOutput3->Add(fhInvMassSumAll5MevBin3);
 415
 416   fhInvMassD0Sum10MevBin4 = new TH1F("fhInvMassD0Sum10MevBin4", "D0 invariant mass Bin4 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
 417   fhInvMassD0Sum10MevBin4->Sumw2(); //Create structure to store sum of squares of weights
 418   fhInvMassD0Sum10MevBin4->SetMinimum(0);
 419   fOutput4->Add(fhInvMassD0Sum10MevBin4);
 420
 421   fhInvMassD0barSum10MevBin4 = new TH1F("fhInvMassD0barSum10MevBin4", "D0bar invariant mass Bin4 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
 422   fhInvMassD0barSum10MevBin4->Sumw2(); //Create structure to store sum of squares of weights
 423   fhInvMassD0barSum10MevBin4->SetMinimum(0);
 424   fOutput4->Add(fhInvMassD0barSum10MevBin4);
 425
 426   fhInvMassSumAll10MevBin4 = new TH1F("fhInvMassSumAll10MevBin4", "D0/D0bar invariant mass Bin4 (good hyps); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
 427   fhInvMassSumAll10MevBin4->Sumw2(); //Create structure to store sum of squares of weights
 428   fhInvMassSumAll10MevBin4->SetMinimum(0);
 429   fOutput4->Add(fhInvMassSumAll10MevBin4);
 430
 431   fhInvMassD0Sum5MevBin4 = new TH1F("fhInvMassD0Sum5MevBin4", "D0 invariant mass Bin4 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
 432   fhInvMassD0Sum5MevBin4->Sumw2(); //Create structure to store sum of squares of weights
 433   fhInvMassD0Sum5MevBin4->SetMinimum(0);
 434   fOutput4->Add(fhInvMassD0Sum5MevBin4);
 435
 436   fhInvMassD0barSum5MevBin4 = new TH1F("fhInvMassD0barSum5MevBin4", "D0bar invariant mass Bin4 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
 437   fhInvMassD0barSum5MevBin4->Sumw2(); //Create structure to store sum of squares of weights
 438   fhInvMassD0barSum5MevBin4->SetMinimum(0);
 439   fOutput4->Add(fhInvMassD0barSum5MevBin4);
 440
 441   fhInvMassSumAll5MevBin4 = new TH1F("fhInvMassSumAll5MevBin4", "D0/D0bar invariant mass Bin4 (good hyps); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
 442   fhInvMassSumAll5MevBin4->Sumw2(); //Create structure to store sum of squares of weights
 443   fhInvMassSumAll5MevBin4->SetMinimum(0);
 444   fOutput4->Add(fhInvMassSumAll5MevBin4);
 445
 446   fhInvMassD0Sum10MevBin5 = new TH1F("fhInvMassD0Sum10MevBin5", "D0 invariant mass Bin5 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
 447   fhInvMassD0Sum10MevBin5->Sumw2(); //Create structure to store sum of squares of weights
 448   fhInvMassD0Sum10MevBin5->SetMinimum(0);
 449   fOutput5->Add(fhInvMassD0Sum10MevBin5);
 450
 451   fhInvMassD0barSum10MevBin5 = new TH1F("fhInvMassD0barSum10MevBin5", "D0bar invariant mass Bin5 (good hyp); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
 452   fhInvMassD0barSum10MevBin5->Sumw2(); //Create structure to store sum of squares of weights
 453   fhInvMassD0barSum10MevBin5->SetMinimum(0);
 454   fOutput5->Add(fhInvMassD0barSum10MevBin5);
 455
 456   fhInvMassSumAll10MevBin5 = new TH1F("fhInvMassSumAll10MevBin5", "D0/D0bar invariant mass Bin5 (good hyps); Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2);
 457   fhInvMassSumAll10MevBin5->Sumw2(); //Create structure to store sum of squares of weights
 458   fhInvMassSumAll10MevBin5->SetMinimum(0);
 459   fOutput5->Add(fhInvMassSumAll10MevBin5);
 460
 461   fhInvMassD0Sum5MevBin5 = new TH1F("fhInvMassD0Sum5MevBin5", "D0 invariant mass Bin5 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
 462   fhInvMassD0Sum5MevBin5->Sumw2(); //Create structure to store sum of squares of weights
 463   fhInvMassD0Sum5MevBin5->SetMinimum(0);
 464   fOutput5->Add(fhInvMassD0Sum5MevBin5);
 465
 466   fhInvMassD0barSum5MevBin5 = new TH1F("fhInvMassD0barSum5MevBin5", "D0bar invariant mass Bin5 (good hyp); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
 467   fhInvMassD0barSum5MevBin5->Sumw2(); //Create structure to store sum of squares of weights
 468   fhInvMassD0barSum5MevBin5->SetMinimum(0);
 469   fOutput5->Add(fhInvMassD0barSum5MevBin5);
 470
 471   fhInvMassSumAll5MevBin5 = new TH1F("fhInvMassSumAll5MevBin5", "D0/D0bar invariant mass Bin5 (good hyps); Inv. mass [GeV]; Entries/5 MeV",120,1.6,2.2);
 472   fhInvMassSumAll5MevBin5->Sumw2(); //Create structure to store sum of squares of weights
 473   fhInvMassSumAll5MevBin5->SetMinimum(0);
 474   fOutput5->Add(fhInvMassSumAll5MevBin5);
 475
 476   fhReflBin1 = new TH2F("fhReflBin1", "Invariant Mass Histogram for reflections; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,8,0.,8.);
 477   fhReflBin1->Sumw2(); //Create structure to store sum of squares of weights
 478   fhReflBin1->SetMinimum(0);
 479   fOutput->Add(fhReflBin1);
 480
 481   fhReflBin2 = new TH2F("fhReflBin2", "Invariant Mass Histogram for reflections; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,8,0.,8.);
 482   fhReflBin2->Sumw2(); //Create structure to store sum of squares of weights
 483   fhReflBin2->SetMinimum(0);
 484   fOutput2->Add(fhReflBin2);
 485
 486   fhReflBin3 = new TH2F("fhReflBin3", "Invariant Mass Histogram for reflections; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,8,0.,8.);
 487   fhReflBin3->Sumw2(); //Create structure to store sum of squares of weights
 488   fhReflBin3->SetMinimum(0);
 489   fOutput3->Add(fhReflBin3);
 490
 491   fhReflBin4 = new TH2F("fhReflBin4", "Invariant Mass Histogram for reflections; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,8,0.,8.);
 492   fhReflBin4->Sumw2(); //Create structure to store sum of squares of weights
 493   fhReflBin4->SetMinimum(0);
 494   fOutput4->Add(fhReflBin4);
 495
 496   fhReflBin5 = new TH2F("fhReflBin5", "Invariant Mass Histogram for reflections; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,8,0.,8.);
 497   fhReflBin5->Sumw2(); //Create structure to store sum of squares of weights
 498   fhReflBin5->SetMinimum(0);
 499   fOutput5->Add(fhReflBin5);
 500
 501   fhReflD0Bin1 = new TH2F("fhReflD0Bin1", "Invariant Mass Histogram for reflections - D0 hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
 502   fhReflD0Bin1->Sumw2(); //Create structure to store sum of squares of weights
 503   fhReflD0Bin1->SetMinimum(0);
 504   fOutput->Add(fhReflD0Bin1);
 505
 506   fhReflD0Bin2 = new TH2F("fhReflD0Bin2", "Invariant Mass Histogram for reflections - D0 hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
 507   fhReflD0Bin2->Sumw2(); //Create structure to store sum of squares of weights
 508   fhReflD0Bin2->SetMinimum(0);
 509   fOutput2->Add(fhReflD0Bin2);
 510
 511   fhReflD0Bin3 = new TH2F("fhReflD0Bin3", "Invariant Mass Histogram for reflections - D0 hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
 512   fhReflD0Bin3->Sumw2(); //Create structure to store sum of squares of weights
 513   fhReflD0Bin3->SetMinimum(0);
 514   fOutput3->Add(fhReflD0Bin3);
 515
 516   fhReflD0Bin4 = new TH2F("fhReflD0Bin4", "Invariant Mass Histogram for reflections - D0 hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
 517   fhReflD0Bin4->Sumw2(); //Create structure to store sum of squares of weights
 518   fhReflD0Bin4->SetMinimum(0);
 519   fOutput4->Add(fhReflD0Bin4);
 520
 521   fhReflD0Bin5 = new TH2F("fhReflD0Bin5", "Invariant Mass Histogram for reflections - D0 hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
 522   fhReflD0Bin5->Sumw2(); //Create structure to store sum of squares of weights
 523   fhReflD0Bin5->SetMinimum(0);
 524   fOutput5->Add(fhReflD0Bin5);
 525
 526   fhReflD0barBin1 = new TH2F("fhReflD0barBin1", "Invariant Mass Histogram for reflections - D0bar hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
 527   fhReflD0barBin1->Sumw2(); //Create structure to store sum of squares of weights
 528   fhReflD0barBin1->SetMinimum(0);
 529   fOutput->Add(fhReflD0barBin1);
 530
 531   fhReflD0barBin2 = new TH2F("fhReflD0barBin2", "Invariant Mass Histogram for reflections - D0bar hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
 532   fhReflD0barBin2->Sumw2(); //Create structure to store sum of squares of weights
 533   fhReflD0barBin2->SetMinimum(0);
 534   fOutput2->Add(fhReflD0barBin2);
 535
 536   fhReflD0barBin3 = new TH2F("fhReflD0barBin3", "Invariant Mass Histogram for reflections - D0bar hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
 537   fhReflD0barBin3->Sumw2(); //Create structure to store sum of squares of weights
 538   fhReflD0barBin3->SetMinimum(0);
 539   fOutput3->Add(fhReflD0barBin3);
 540
 541   fhReflD0barBin4 = new TH2F("fhReflD0barBin4", "Invariant Mass Histogram for reflections - D0bar hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
 542   fhReflD0barBin4->Sumw2(); //Create structure to store sum of squares of weights
 543   fhReflD0barBin4->SetMinimum(0);
 544   fOutput4->Add(fhReflD0barBin4);
 545
 546   fhReflD0barBin5 = new TH2F("fhReflD0barBin5", "Invariant Mass Histogram for reflections - D0bar hyp.; Inv. mass [GeV]; Entries/10 MeV",60,1.6,2.2,7,0.,7.);
 547   fhReflD0barBin5->Sumw2(); //Create structure to store sum of squares of weights
 548   fhReflD0barBin5->SetMinimum(0);
 549   fOutput5->Add(fhReflD0barBin5);
 550
 551   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.);
 552   fScatterP4PID->SetMinimum(0);
 553   fOutput->Add(fScatterP4PID);
 554
 555   fPtVsY = new TH2F("fPtVsY", "Pt vs Y PPR Sel. Candidates; Pt [GeV/c]; Y",250,0.,25.,300,-3.,3.);
 556   fPtVsY->SetMinimum(0);
 557   fOutputC->Add(fPtVsY);
 558
 559   fPtVsYAll = new TH2F("fPtVsYAll", "Pt vs Y All Candidates; Pt [GeV/c]; Y",250,0.,25.,300,-3.,3.);
 560   fPtVsYAll->SetMinimum(0);
 561   fOutputC->Add(fPtVsYAll);
 562
 563   fEventCounter = new TH1F("fEventCounter", "N° of total events; NA; Events",1,0.,1.);
 564   fEventCounter->SetMinimum(0);
 565   fOutputC->Add(fEventCounter);
 566
 567   fCutDCA = new TH1F("fCutDCA", "DCA of candidate (couple); DCA [cm]; Entries/micron",500,0.,0.05);
 568   fCutDCA->SetMinimum(0);
 569   fOutputC->Add(fCutDCA);
 570
 571   fCutDCA3 = new TH1F("fCutDCA3", "DCA of candidate (trips); DCA [cm]; Entries/micron",500,0.,0.05);
 572   fCutDCA3->SetMinimum(0);
 573   fOutputC->Add(fCutDCA3);
 574
 575   fCutDCA2 = new TH1F("fCutDCA2", "DCA of candidate (quads1); DCA [cm]; Entries/micron",500,0.,0.05);
 576   fCutDCA2->SetMinimum(0);
 577   fOutputC->Add(fCutDCA2);
 578
 579   fCutDCA5 = new TH1F("fCutDCA5", "DCA of candidate (quads2); DCA [cm]; Entries/micron",500,0.,0.05);
 580   fCutDCA5->SetMinimum(0);
 581   fOutputC->Add(fCutDCA5);
 582
 583   fCutVertexDist2 = new TH1F("fCutVertexDist2", "Distance Vtx doubl.-Primary Vtx; Distance [cm]; Entries/15 micron",500,0.,0.75);
 584   fCutVertexDist2->SetMinimum(0);
 585   fOutputC->Add(fCutVertexDist2);
 586
 587   fCutVertexDist3 = new TH1F("fCutVertexDist3", "Distance Vtx trips-Primary Vtx; Distance [cm]; Entries/10 micron",500,0.,0.5);
 588   fCutVertexDist3->SetMinimum(0);
 589   fOutputC->Add(fCutVertexDist3);
 590
 591   fCutVertexDist4 = new TH1F("fCutVertexDist4", "Distance Vtx quads-Primary Vtx; Distance [cm]; Entries/5 micron",500,0.,0.25);
 592   fCutVertexDist4->SetMinimum(0);
 593   fOutputC->Add(fCutVertexDist4);
 594
 595   fCutCosinePoint = new TH1F("fCutCosinePoint", "Cosine of angle of pointing; Cos(Thetapt.); Entries/10^(-3)",250,0.75,1.);
 596   fCutCosinePoint->SetMinimum(0);
 597   fOutputC->Add(fCutCosinePoint);
 598
 599   fCutPt = new TH1F("fCutPt", "Pt of candidate D0; Pt [GeV/c]; Entries/5 MeV",3000,0.,15.);
 600   fCutPt->SetMinimum(0);
 601   fOutputC->Add(fCutPt);
 602
 603   fCutY = new TH1F("fCutY", "Y of candidate D0; Pt [GeV/c]; Entries/5 MeV",900,-9.,9.);
 604   fCutY->SetMinimum(0);
 605   fOutputC->Add(fCutY);
 606
 607   fPIDSel = new TH1F("fPIDSel", "Ratio of D0 selected by PID for Correction",3,0.,3.);
 608   fPIDSel->SetMinimum(0);
 609   fPIDSel->GetXaxis()->SetBinLabel(1,"D0allhyp All");
 610   fPIDSel->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
 611   fPIDSel->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
 612
 613   fPIDSelBin1 = new TH1F("fPIDSelBin1", "Ratio of D0 selected by PID for Correction",3,0.,3.);
 614   fPIDSelBin1->SetMinimum(0);
 615   fPIDSelBin1->GetXaxis()->SetBinLabel(1,"D0allhyp All");
 616   fPIDSelBin1->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
 617   fPIDSelBin1->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
 618
 619   fPIDSelBin2 = new TH1F("fPIDSelBin1", "Ratio of D0 selected by PID for Correction",3,0.,3.);
 620   fPIDSelBin2->SetMinimum(0);
 621   fPIDSelBin2->GetXaxis()->SetBinLabel(1,"D0allhyp All");
 622   fPIDSelBin2->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
 623   fPIDSelBin2->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
 624
 625   fPIDSelBin3 = new TH1F("fPIDSelBin1", "Ratio of D0 selected by PID for Correction",3,0.,3.);
 626   fPIDSelBin3->SetMinimum(0);
 627   fPIDSelBin3->GetXaxis()->SetBinLabel(1,"D0allhyp All");
 628   fPIDSelBin3->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
 629   fPIDSelBin3->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
 630
 631   fPIDSelBin4 = new TH1F("fPIDSelBin1", "Ratio of D0 selected by PID for Correction",3,0.,3.);
 632   fPIDSelBin4->SetMinimum(0);
 633   fPIDSelBin4->GetXaxis()->SetBinLabel(1,"D0allhyp All");
 634   fPIDSelBin4->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
 635   fPIDSelBin4->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
 636
 637   fPIDSelBin5 = new TH1F("fPIDSelBin1", "Ratio of D0 selected by PID for Correction",3,0.,3.);
 638   fPIDSelBin5->SetMinimum(0);
 639   fPIDSelBin5->GetXaxis()->SetBinLabel(1,"D0allhyp All");
 640   fPIDSelBin5->GetXaxis()->SetBinLabel(2,"D0allhyp PID");
 641   fPIDSelBin5->GetXaxis()->SetBinLabel(3,"D0allhyp PID (hypok)");
 642
 643   fMultipleHyps = new TH2F("fMultipleHyps", "N. of hyp. accepted for each candidate (accounted N. times)",8,0.,8.,5,0.,5.);
 644   fMultipleHyps->SetMinimum(0);
 645   fMultipleHyps->GetXaxis()->SetBinLabel(1,"1 (PPR)");
 646   fMultipleHyps->GetXaxis()->SetBinLabel(2,"2 (PPR)");
 647   fMultipleHyps->GetXaxis()->SetBinLabel(3,"3 (PPR)");
 648   fMultipleHyps->GetXaxis()->SetBinLabel(4,"4 (PPR)");
 649   fMultipleHyps->GetXaxis()->SetBinLabel(5,"1 (PPR+PID)");
 650   fMultipleHyps->GetXaxis()->SetBinLabel(6,"2 (PPR+PID)");
 651   fMultipleHyps->GetXaxis()->SetBinLabel(7,"3 (PPR+PID)");
 652   fMultipleHyps->GetXaxis()->SetBinLabel(8,"4 (PPR+PID)");
 653   fMultipleHyps->GetYaxis()->SetBinLabel(1,"PtBin 1");
 654   fMultipleHyps->GetYaxis()->SetBinLabel(2,"PtBin 2");
 655   fMultipleHyps->GetYaxis()->SetBinLabel(3,"PtBin 3");
 656   fMultipleHyps->GetYaxis()->SetBinLabel(4,"PtBin 4");
 657   fMultipleHyps->GetYaxis()->SetBinLabel(5,"PtBin 5");
 658
 659   fMultipleHypsType = new TH2F("fMultipleHypsType", "Type of hyp. accepted for each candidate",8,0.,8.,5,0.,5.);
 660   fMultipleHypsType->SetMinimum(0);
 661   fMultipleHypsType->GetXaxis()->SetBinLabel(1,"D0");
 662   fMultipleHypsType->GetXaxis()->SetBinLabel(2,"D0bar");
 663   fMultipleHypsType->GetXaxis()->SetBinLabel(3,"2D0");
 664   fMultipleHypsType->GetXaxis()->SetBinLabel(4,"2D0bar");
 665   fMultipleHypsType->GetXaxis()->SetBinLabel(5,"D0+D0bar");
 666   fMultipleHypsType->GetXaxis()->SetBinLabel(6,"2D0+D0bar");
 667   fMultipleHypsType->GetXaxis()->SetBinLabel(7,"D0+2D0bar");
 668   fMultipleHypsType->GetXaxis()->SetBinLabel(8,"2D0+2D0bar");
 669   fMultipleHypsType->GetYaxis()->SetBinLabel(1,"PtBin 1");
 670   fMultipleHypsType->GetYaxis()->SetBinLabel(2,"PtBin 2");
 671   fMultipleHypsType->GetYaxis()->SetBinLabel(3,"PtBin 3");
 672   fMultipleHypsType->GetYaxis()->SetBinLabel(4,"PtBin 4");
 673   fMultipleHypsType->GetYaxis()->SetBinLabel(5,"PtBin 5");
 674
 675   fOutputC->Add(fMultipleHyps);
 676   fOutputC->Add(fMultipleHypsType);
 677
 678   fOutputC->Add(fPIDSel);
 679   fOutput->Add(fPIDSelBin1);
 680   fOutput2->Add(fPIDSelBin2);
 681   fOutput3->Add(fPIDSelBin3);
 682   fOutput4->Add(fPIDSelBin4);
 683   fOutput5->Add(fPIDSelBin5);
 684
 685   fPtSel = new TH1F("fPtSel", "Pt of candidates accepted; Pt [GeV/c]; Entries/10 MeV",2000,0.,20.);
 686   fPtSel->SetMinimum(0);
 687   fOutputC->Add(fPtSel);
 688
 689   return;
 690 }
 691
 692 //________________________________________________________________________
 693 void AliAnalysisTaskSESelectHF4Prong::UserExec(Option_t */*option*/)
 694 {
 695   // Execute analysis for current event:
 696   // heavy flavor candidates selection and histograms
 697
 698   AliAODEvent *aodIn = dynamic_cast<AliAODEvent*> (InputEvent());
 699
 700   TClonesArray *inputArrayCharm4Prong = 0;
 701
 702   if(!aodIn && AODEvent() && IsStandardAOD()) {
 703     // In case there is an AOD handler writing a standard AOD, use the AOD
 704     // event in memory rather than the input (ESD) event.
 705     aodIn = dynamic_cast<AliAODEvent*> (AODEvent());
 706     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 707     // have to taken from the AOD event hold by the AliAODExtension
 708     AliAODHandler* aodHandler = (AliAODHandler*)
 709       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 710     if(aodHandler->GetExtensions()) {
 711       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 712       AliAODEvent *aodFromExt = ext->GetAOD();
 713       // load D0 candidates
 714       inputArrayCharm4Prong=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
 715     }
 716   } else {
 717     // load D0 candidates
 718     inputArrayCharm4Prong=(TClonesArray*)aodIn->GetList()->FindObject("Charm4Prong");
 719   }
 720
 721   if(!inputArrayCharm4Prong) {
 722     printf("AliAnalysisTaskSESelectHF4Prong::UserExec: D0to3Kpi branch not found!\n");
 723     return;
 724   }
 725
 726   //print event info
 727 //  aodIn->GetHeader()->Print();
 728
 729   //Event counter ++
 730   fEventCounter->Fill(0);
 731
 732   // fix for temporary bug in ESDfilter
 733   // the AODs with null vertex pointer didn't pass the PhysSel
 734   if(!aodIn->GetPrimaryVertex()) return;
 735
 736   // primary vertex
 737   AliAODVertex *vtx1 = (AliAODVertex*)aodIn->GetPrimaryVertex();
 738 //  vtx1->Print();
 739
 740   // make trkIDtoEntry register (temporary)
 741   Int_t trkIDtoEntry[100000];
 742   for(Int_t it=0;it<aodIn->GetNumberOfTracks();it++) {
 743     AliAODTrack *track = aodIn->GetTrack(it);
 744     trkIDtoEntry[track->GetID()]=it;
 745   }
 746
 747   Int_t iOutVerticesHF=0,iOutCharm4Prong=0;
 748   fVerticesHFTClArr->Delete();
 749   iOutVerticesHF = fVerticesHFTClArr->GetEntriesFast();
 750   TClonesArray &verticesHFRef = *fVerticesHFTClArr;
 751   fCharm4ProngTClArr->Delete();
 752   iOutCharm4Prong = fCharm4ProngTClArr->GetEntriesFast();
 753   TClonesArray &aodCharm4ProngRef = *fCharm4ProngTClArr;
 754
 755   // loop over D0->K3pi candidates
 756   Int_t nInCharm4Prong = inputArrayCharm4Prong->GetEntriesFast();
 757   printf("Number of D0->K3pi: %d\n",nInCharm4Prong);
 758
 759   for (Int_t iCharm4Prong = 0; iCharm4Prong < nInCharm4Prong; iCharm4Prong++) {
 760     AliAODRecoDecayHF4Prong *dIn = (AliAODRecoDecayHF4Prong*)inputArrayCharm4Prong->UncheckedAt(iCharm4Prong);
 761     Bool_t unsetvtx=kFALSE;
 762
 763     if(!dIn->GetOwnPrimaryVtx()) {
 764       dIn->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
 765       unsetvtx=kTRUE;
 766     }
 767
 768    //fill histos of cuts
 769    Double_t dca = dIn->GetDCA();
 770    Double_t dist2 = dIn->GetDist12toPrim();
 771    Double_t dist3 = dIn->GetDist3toPrim();
 772    Double_t dist4 = dIn->GetDist4toPrim();
 773    Double_t cosine = dIn->CosPointingAngle();
 774    Double_t ptPart = dIn->Pt();
 775    Double_t yPart = dIn->YD0();
 776    Double_t dcatrip = dIn->GetDCA(3);
 777    Double_t dcaquad1 = dIn->GetDCA(2);
 778    Double_t dcaquad2 = dIn->GetDCA(5);
 779
 780    fCutDCA->Fill(dca);
 781    fCutDCA3->Fill(dcatrip);
 782    fCutDCA2->Fill(dcaquad1);
 783    fCutDCA5->Fill(dcaquad2);
 784    fCutVertexDist2->Fill(dist2);
 785    fCutVertexDist3->Fill(dist3);
 786    fCutVertexDist4->Fill(dist4);
 787    fCutCosinePoint->Fill(cosine);
 788    fCutPt->Fill(ptPart);
 789    fCutY->Fill(yPart);
 790    fPtVsYAll->Fill(ptPart,yPart);
 791
 792    //flags initialization
 793    fSelected = fCuts->IsSelected(dIn,AliRDHFCuts::kCandidate);
 794    Int_t selD01 = fCuts->D01Selected(dIn,AliRDHFCuts::kCandidate);
 795    Int_t selD02 = fCuts->D02Selected(dIn,AliRDHFCuts::kCandidate);
 796    Int_t selD0bar1 = fCuts->D0bar1Selected(dIn,AliRDHFCuts::kCandidate);
 797    Int_t selD0bar2 = fCuts->D0bar2Selected(dIn,AliRDHFCuts::kCandidate);
 798    Int_t flagAccLim = 1;
 799
 800    //Limited Acceptance
 801    if(ptPart > 5.) {
 802      if (TMath::Abs(yPart) > 0.8) flagAccLim = 0;
 803    }
 804    else {
 805      Double_t maxFiducialY = -0.2/15*ptPart*ptPart+1.9/15*ptPart+0.5;
 806      Double_t minFiducialY = 0.2/15*ptPart*ptPart-1.9/15*ptPart-0.5;
 807      if (yPart < minFiducialY || yPart > maxFiducialY) flagAccLim = 0;;
 808    }
 809
 810    //number of CANDIDATES (regardless of hypotheses) passing PPR
 811    if(fSelected==1||fSelected==2||fSelected==3) {
 812       fPtVsY->Fill(ptPart,yPart);
 813       fPIDSel->Fill(0);
 814         if (ptPart >= fPtBinH[0] && ptPart < fPtBinH[1])      fPIDSelBin1->Fill(0);
 815         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fPIDSelBin2->Fill(0);
 816         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fPIDSelBin3->Fill(0);
 817         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fPIDSelBin4->Fill(0);
 818         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fPIDSelBin5->Fill(0);
 819         }
 820
 821    PostData(6,fOutputC);
 822
 823    //selection
 824    if((fSelected==1||fSelected==2||fSelected==3) && flagAccLim == 1) {
 825       // get daughter AOD tracks
 826       AliAODTrack *trk0 = (AliAODTrack*)dIn->GetDaughter(0);
 827       AliAODTrack *trk1 = (AliAODTrack*)dIn->GetDaughter(1);
 828       AliAODTrack *trk2 = (AliAODTrack*)dIn->GetDaughter(2);
 829       AliAODTrack *trk3 = (AliAODTrack*)dIn->GetDaughter(3);
 830       if(!trk0 || !trk1 || !trk2 || !trk3) {
 831         trk0=aodIn->GetTrack(trkIDtoEntry[dIn->GetProngID(0)]);
 832         trk1=aodIn->GetTrack(trkIDtoEntry[dIn->GetProngID(1)]);
 833         trk2=aodIn->GetTrack(trkIDtoEntry[dIn->GetProngID(2)]);
 834         trk3=aodIn->GetTrack(trkIDtoEntry[dIn->GetProngID(3)]);
 835       }
 836       printf("pt of positive track #1: %f\n",trk0->Pt());
 837       printf("pt of negative track #1: %f\n",trk1->Pt());
 838       printf("pt of positive track #2: %f\n",trk2->Pt());
 839       printf("pt of negative track #2: %f\n",trk3->Pt());
 840
 841       dIn->InvMassD0(fmassD0);
 842       dIn->InvMassD0bar(fmassD0bar);
 843
 844       //fill histos (combining selection form cuts & PID (with rho information))
 845       Double_t binPt = -1;
 846       Int_t hypD01 = 0, hypD02 = 0, hypD0bar1 = 0, hypD0bar2 = 0;
 847       Int_t pid1 = 0, pid2 = 0, pidbar1 = 0, pidbar2 = 0;
 848       Int_t pidSelection = fCuts->IsSelectedFromPID(dIn, &pid1, &pid2, &pidbar1, &pidbar2);
 849
 850    //number of CANDIDATES (regardless of hypotheses) passing PPR + PID - PAY ATTENTION: hypoth. for PID and PPR may not be the same!
 851    if (pidSelection > 0) {
 852         fPIDSel->Fill(1);
 853         if (ptPart >= fPtBinH[0] && ptPart < fPtBinH[1])      {fPIDSelBin1->Fill(1); binPt = 0;}
 854         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) {fPIDSelBin2->Fill(1); binPt = 1;}
 855         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) {fPIDSelBin3->Fill(1); binPt = 2;}
 856         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) {fPIDSelBin4->Fill(1); binPt = 3;}
 857         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) {fPIDSelBin5->Fill(1); binPt = 4;}
 858       }
 859
 860            //number of hypoteses accepted per candidate after PPR
 861            if(selD01+selD02+selD0bar1+selD0bar2 == 1) fMultipleHyps->Fill((Double_t)0,binPt);
 862            if(selD01+selD02+selD0bar1+selD0bar2 == 2) fMultipleHyps->Fill((Double_t)1,binPt);
 863            if(selD01+selD02+selD0bar1+selD0bar2 == 3) fMultipleHyps->Fill((Double_t)2,binPt);
 864            if(selD01+selD02+selD0bar1+selD0bar2 == 4) fMultipleHyps->Fill((Double_t)3,binPt);
 865
 866    //combine PPD + PID cuts
 867    if (selD01 == 1 && pid1==1) hypD01 = 1;
 868    if (selD02 == 1 && pid2==1) hypD02 = 1;
 869    if (selD0bar1 == 1 && pidbar1==1) hypD0bar1 = 1;
 870    if (selD0bar2 == 1 && pidbar2==1) hypD0bar2 = 1;
 871
 872    //number of CANDIDATES (regardless of hypotheses) passing PPR + PID - PAY ATTENTION: hypoth. for PID and PPR must match (at least one)!
 873    if (hypD01 == 1 || hypD02 == 1 || hypD0bar1 == 1 || hypD0bar2 == 1) {
 874         fPIDSel->Fill(2);
 875         if (ptPart >= fPtBinH[0] && ptPart < fPtBinH[1])      {fPIDSelBin1->Fill(2); binPt = 0;}
 876         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) {fPIDSelBin2->Fill(2); binPt = 1;}
 877         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) {fPIDSelBin3->Fill(2); binPt = 2;}
 878         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) {fPIDSelBin4->Fill(2); binPt = 3;}
 879         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) {fPIDSelBin5->Fill(2); binPt = 4;}
 880       }
 881
 882            //number of hypoteses accepted per candidate after PPR and PID
 883            if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 1) fMultipleHyps->Fill((Double_t)4,binPt);
 884            if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 2) fMultipleHyps->Fill((Double_t)5,binPt);
 885            if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 3) fMultipleHyps->Fill((Double_t)6,binPt);
 886            if(hypD01+hypD02+hypD0bar1+hypD0bar2 == 4) fMultipleHyps->Fill((Double_t)7,binPt);
 887
 888            //type of hypoteses accepted per candidate after PPR and PID
 889            if(hypD01+hypD02 == 1 && hypD0bar1+hypD0bar2 == 0) fMultipleHypsType->Fill((Double_t)0,binPt);
 890            if(hypD01+hypD02 == 0 && hypD0bar1+hypD0bar2 == 1) fMultipleHypsType->Fill((Double_t)1,binPt);
 891            if(hypD01+hypD02 == 2 && hypD0bar1+hypD0bar2 == 0) fMultipleHypsType->Fill((Double_t)2,binPt);
 892            if(hypD01+hypD02 == 0 && hypD0bar1+hypD0bar2 == 2) fMultipleHypsType->Fill((Double_t)3,binPt);
 893            if(hypD01+hypD02 == 1 && hypD0bar1+hypD0bar2 == 1) fMultipleHypsType->Fill((Double_t)4,binPt);
 894            if(hypD01+hypD02 == 2 && hypD0bar1+hypD0bar2 == 1) fMultipleHypsType->Fill((Double_t)5,binPt);
 895            if(hypD01+hypD02 == 1 && hypD0bar1+hypD0bar2 == 2) fMultipleHypsType->Fill((Double_t)6,binPt);
 896            if(hypD01+hypD02 == 2 && hypD0bar1+hypD0bar2 == 2) fMultipleHypsType->Fill((Double_t)7,binPt);
 897
 898    //Call function for reflection analysis
 899    if (hypD01+hypD02+hypD0bar1+hypD0bar2 > 0) AnalysisReflection(aodIn, dIn, hypD01, hypD02, hypD0bar1, hypD0bar2);
 900
 901    //All histos are filled if Pt of candidate is greater than minimum of first bin (in this way: bin1+bin2+...binN = whole)
 902    if (ptPart > fPtBinH[0]) {
 903
 904       // D01 hyp.
 905       if(hypD01==1) {
 906         fPtSel->Fill(ptPart);
 907         fScatterP4PID->Fill(trk1->Pt(),trk3->Pt());
 908         if (ptPart < fPtBinH[1])                          {fhInvMassD0Sum10MevBin1->Fill(fmassD0[0]);
 909                                                           fhInvMassD0Sum5MevBin1->Fill(fmassD0[0]);
 910                                                           fhInvMassSumAll10MevBin1->Fill(fmassD0[0]);
 911                                                           fhInvMassSumAll5MevBin1->Fill(fmassD0[0]);}
 912         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) {fhInvMassD0Sum10MevBin2->Fill(fmassD0[0]);
 913                                                               fhInvMassD0Sum5MevBin2->Fill(fmassD0[0]);
 914                                                               fhInvMassSumAll10MevBin2->Fill(fmassD0[0]);
 915                                                               fhInvMassSumAll5MevBin2->Fill(fmassD0[0]);}
 916         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) {fhInvMassD0Sum10MevBin3->Fill(fmassD0[0]);
 917                                                               fhInvMassD0Sum5MevBin3->Fill(fmassD0[0]);
 918                                                               fhInvMassSumAll10MevBin3->Fill(fmassD0[0]);
 919                                                               fhInvMassSumAll5MevBin3->Fill(fmassD0[0]);}
 920         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) {fhInvMassD0Sum10MevBin4->Fill(fmassD0[0]);
 921                                                               fhInvMassD0Sum5MevBin4->Fill(fmassD0[0]);
 922                                                               fhInvMassSumAll10MevBin4->Fill(fmassD0[0]);
 923                                                               fhInvMassSumAll5MevBin4->Fill(fmassD0[0]);}
 924         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5])   {fhInvMassD0Sum10MevBin5->Fill(fmassD0[0]);
 925                                                                 fhInvMassD0Sum5MevBin5->Fill(fmassD0[0]);
 926                                                                 fhInvMassSumAll10MevBin5->Fill(fmassD0[0]);
 927                                                                 fhInvMassSumAll5MevBin5->Fill(fmassD0[0]);}
 928      }
 929       // D02 hyp.
 930       if(hypD02==1) {
 931         fPtSel->Fill(ptPart);
 932         fScatterP4PID->Fill(trk3->Pt(),trk1->Pt());
 933         if (ptPart < fPtBinH[1])                          {fhInvMassD0Sum10MevBin1->Fill(fmassD0[1]);
 934                                                           fhInvMassD0Sum5MevBin1->Fill(fmassD0[1]);
 935                                                           fhInvMassSumAll10MevBin1->Fill(fmassD0[1]);
 936                                                           fhInvMassSumAll5MevBin1->Fill(fmassD0[1]);}
 937         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) {fhInvMassD0Sum10MevBin2->Fill(fmassD0[1]);
 938                                                               fhInvMassD0Sum5MevBin2->Fill(fmassD0[1]);
 939                                                               fhInvMassSumAll10MevBin2->Fill(fmassD0[1]);
 940                                                               fhInvMassSumAll5MevBin2->Fill(fmassD0[1]);}
 941         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) {fhInvMassD0Sum10MevBin3->Fill(fmassD0[1]);
 942                                                               fhInvMassD0Sum5MevBin3->Fill(fmassD0[1]);
 943                                                               fhInvMassSumAll10MevBin3->Fill(fmassD0[1]);
 944                                                               fhInvMassSumAll5MevBin3->Fill(fmassD0[1]);}
 945         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) {fhInvMassD0Sum10MevBin4->Fill(fmassD0[1]);
 946                                                               fhInvMassD0Sum5MevBin4->Fill(fmassD0[1]);
 947                                                               fhInvMassSumAll10MevBin4->Fill(fmassD0[1]);
 948                                                               fhInvMassSumAll5MevBin4->Fill(fmassD0[1]);}
 949         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5])   {fhInvMassD0Sum10MevBin5->Fill(fmassD0[1]);
 950                                                                  fhInvMassD0Sum5MevBin5->Fill(fmassD0[1]);
 951                                                                  fhInvMassSumAll10MevBin5->Fill(fmassD0[1]);
 952                                                                  fhInvMassSumAll5MevBin5->Fill(fmassD0[1]);}
 953      }
 954       // D0bar1 hyp.
 955       if(hypD0bar1==1) {
 956         fPtSel->Fill(ptPart);
 957         fScatterP4PID->Fill(trk0->Pt(),trk2->Pt());
 958         if (ptPart < fPtBinH[1])                          {fhInvMassD0barSum10MevBin1->Fill(fmassD0bar[0]);
 959                                                           fhInvMassD0barSum5MevBin1->Fill(fmassD0bar[0]);
 960                                                           fhInvMassSumAll10MevBin1->Fill(fmassD0bar[0]);
 961                                                           fhInvMassSumAll5MevBin1->Fill(fmassD0bar[0]);}
 962         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) {fhInvMassD0barSum10MevBin2->Fill(fmassD0bar[0]);
 963                                                               fhInvMassD0barSum5MevBin2->Fill(fmassD0bar[0]);
 964                                                               fhInvMassSumAll10MevBin2->Fill(fmassD0bar[0]);
 965                                                               fhInvMassSumAll5MevBin2->Fill(fmassD0bar[0]);}
 966         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) {fhInvMassD0barSum10MevBin3->Fill(fmassD0bar[0]);
 967                                                               fhInvMassD0barSum5MevBin3->Fill(fmassD0bar[0]);
 968                                                               fhInvMassSumAll10MevBin3->Fill(fmassD0bar[0]);
 969                                                               fhInvMassSumAll5MevBin3->Fill(fmassD0bar[0]);}
 970         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) {fhInvMassD0barSum10MevBin4->Fill(fmassD0bar[0]);
 971                                                               fhInvMassD0barSum5MevBin4->Fill(fmassD0bar[0]);
 972                                                               fhInvMassSumAll10MevBin4->Fill(fmassD0bar[0]);
 973                                                               fhInvMassSumAll5MevBin4->Fill(fmassD0bar[0]);}
 974         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5])   {fhInvMassD0barSum10MevBin5->Fill(fmassD0bar[0]);
 975                                                                 fhInvMassD0barSum5MevBin5->Fill(fmassD0bar[0]);
 976                                                                 fhInvMassSumAll10MevBin5->Fill(fmassD0bar[0]);
 977                                                                 fhInvMassSumAll5MevBin5->Fill(fmassD0bar[0]);}
 978      }
 979       // D0bar2 hyp.
 980       if(hypD0bar2==1) {
 981         fPtSel->Fill(ptPart);
 982         fScatterP4PID->Fill(trk2->Pt(),trk0->Pt());
 983         if (ptPart < fPtBinH[1])                          {fhInvMassD0barSum10MevBin1->Fill(fmassD0bar[1]);
 984                                                           fhInvMassD0barSum5MevBin1->Fill(fmassD0bar[1]);
 985                                                           fhInvMassSumAll10MevBin1->Fill(fmassD0bar[1]);
 986                                                           fhInvMassSumAll5MevBin1->Fill(fmassD0bar[1]);}
 987         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) {fhInvMassD0barSum10MevBin2->Fill(fmassD0bar[1]);
 988                                                               fhInvMassD0barSum5MevBin2->Fill(fmassD0bar[1]);
 989                                                               fhInvMassSumAll10MevBin2->Fill(fmassD0bar[1]);
 990                                                               fhInvMassSumAll5MevBin2->Fill(fmassD0bar[1]);}
 991         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) {fhInvMassD0barSum10MevBin3->Fill(fmassD0bar[1]);
 992                                                               fhInvMassD0barSum5MevBin3->Fill(fmassD0bar[1]);
 993                                                               fhInvMassSumAll10MevBin3->Fill(fmassD0bar[1]);
 994                                                               fhInvMassSumAll5MevBin3->Fill(fmassD0bar[1]);}
 995         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) {fhInvMassD0barSum10MevBin4->Fill(fmassD0bar[1]);
 996                                                               fhInvMassD0barSum5MevBin4->Fill(fmassD0bar[1]);
 997                                                               fhInvMassSumAll10MevBin4->Fill(fmassD0bar[1]);
 998                                                               fhInvMassSumAll5MevBin4->Fill(fmassD0bar[1]);}
 999         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5])   {fhInvMassD0barSum10MevBin5->Fill(fmassD0bar[1]);
1000                                                                 fhInvMassD0barSum5MevBin5->Fill(fmassD0bar[1]);
1001                                                                 fhInvMassSumAll10MevBin5->Fill(fmassD0bar[1]);
1002                                                                 fhInvMassSumAll5MevBin5->Fill(fmassD0bar[1]);}
1003      }
1004    }
1005       PostData(1,fOutput);
1006       PostData(2,fOutput2);
1007       PostData(3,fOutput3);
1008       PostData(4,fOutput4);
1009       PostData(5,fOutput5);
1010
1011       // HERE ONE COULD RECALCULATE THE VERTEX USING THE KF PACKAGE
1012
1013       // clone candidate for output AOD
1014       if(hypD01||hypD02||hypD0bar1||hypD0bar2) {
1015       AliAODVertex *v = new(verticesHFRef[iOutVerticesHF++])
1016         AliAODVertex(*(dIn->GetSecondaryVtx()));
1017       AliAODRecoDecayHF4Prong *dOut=new(aodCharm4ProngRef[iOutCharm4Prong++])
1018         AliAODRecoDecayHF4Prong(*dIn);
1019       dOut->SetSecondaryVtx(v);
1020       dOut->SetOwnPrimaryVtx((AliAODVertex*)((dIn->GetOwnPrimaryVtx())->Clone()));
1021       v->SetParent(dOut);
1022       }
1023     }
1024     if(unsetvtx) dIn->UnsetOwnPrimaryVtx();
1025   } // end loop on D0->K3pi
1026
1027   printf("Number of selected D0->K3pi: %d\n",iOutCharm4Prong);
1028
1029   return;
1030 }
1031
1032 //________________________________________________________________________
1033 void AliAnalysisTaskSESelectHF4Prong::AnalysisReflection(AliAODEvent* aodIn, AliAODRecoDecayHF4Prong* d, Int_t hypD01, Int_t hypD02, Int_t hypD0bar1, Int_t hypD0bar2)
1034 {
1035   /*
1036   ---STUDY OF REFLECTIONS ON CANDIDATE IN ANALYSIS---
1037   Layers of TH2F fhRefl:
1038   0 = all candidates, single hyps selected only;
1039   1 = all candidates, multi hyps selected only;
1040   2 = true D0/D0bar only, all hypotheses selected;
1041   3 = true D0/D0bar only, single hyps selected only; -> I assume that the unique selected hypotesis is the right one!
1042   4 = true D0/D0bar only, multi hyps selected only;
1043   5 = true D0/D0bar only, multi TRUE hyps selected only;
1044   6 = true D0/D0bar only, multi FAKE hyps selected only;
1045   7 = false D0/D0bar only (background)
1046   Layers of TH2F fhReflD0 (idem for D0bar, inverting particles):
1047   0 = true D0 only, true hypothesis for both single and multi hyps cases;
1048   1 = true D0 only, true hypothesis for single hyps case only; -> I assume that the unique selected hypotesis is the right one!
1049   2 = true D0 only, true hypothesis for multi hyps case only;
1050   3 = true D0 only, other D0 wrong hypothesis (multi case, obviously);
1051   4 = true D0 only, D0bar1 wrong hypothesis (multi case, obviously);
1052   5 = true D0 only, D0bar2 wrong hypothesis (multi case, obviously);
1053   6 = true D0 only, D0bar1 + D0bar2 wrong hypotheses
1054   */
1055
1056   Int_t flagMult = hypD01+hypD02+hypD0bar1+hypD0bar2; //single or multi hyps selected
1057   Int_t flagLayer1 = -1, flagLayer2 = -1, flagLayer3 = -1, flagLayer4 = -1; //to select layers in fhRefl
1058
1059   if (flagMult == 1) {flagLayer1 = 0; flagLayer2 = 0; flagLayer3 = 0; flagLayer4 = 0;}
1060   else if (flagMult == 2 || flagMult == 3 || flagMult == 4) {flagLayer1 = 1; flagLayer2 = 1; flagLayer3 = 1; flagLayer4 = 1;}
1061
1062   FillReflHistos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, flagLayer1, flagLayer2, flagLayer3, flagLayer4);  //Fill layers 0 and 1
1063
1064   if (fMCTruth==0) return;
1065
1066   //start of MC Truth phase
1067   TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodIn->FindListObject(AliAODMCParticle::StdBranchName()));
1068   if (!mcArray) {
1069     AliError("Could not find Monte-Carlo in AOD");
1070     return;}
1071
1072   Int_t pdgCand = 421;
1073   Int_t pdgDgCharm4Prong[4]={321,211,211,211}; //pdg of daughters
1074
1075   Int_t mcLabel = d->MatchToMC(pdgCand,mcArray,4,pdgDgCharm4Prong); //selection of true or false candidate (regardless of hypothesis) through MCtruth
1076   printf("MatchToMC = %d\n",mcLabel);
1077
1078   if (mcLabel==-1) { //fill layer 7 (background)
1079   FillReflHistos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, 7, 7, 7, 7);
1080   return;}
1081
1082   FillReflHistos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, 2, 2, 2, 2); //fill layer 2 (true D0/D0bar)
1083
1084   Int_t truthHyp = StudyMCTruth(mcArray, d); //calls function which studies which hypothesis is true for candidate
1085
1086   if (flagMult == 1) { //fill layer 3 (true D0/D0bar - single hyps only)
1087     FillReflHistos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, 3, 3, 3, 3);
1088     switch(truthHyp) { //fills fhReflD0 and fhReflD0bar with layers containing all true D0 or D0bars, for single and all hypotheses (layers 1, 0)
1089       case(1): FillReflD0Histos(d, hypD01, 0, 0, 0, 0, 0, 0, 0);
1090                FillReflD0Histos(d, hypD01, 0, 0, 0, 1, 1, 1, 1); //here I discard the cases in which only a wrong hyp is selected (very very few)
1091                break;
1092       case(2): FillReflD0Histos(d, 0, hypD02, 0, 0, 0, 0, 0, 0);
1093                FillReflD0Histos(d, 0, hypD02, 0, 0, 1, 1, 1, 1);
1094                break;
1095       case(3): FillReflD0barHistos(d, 0, 0, hypD0bar1, 0, 0, 0, 0, 0);
1096                FillReflD0barHistos(d, 0, 0, hypD0bar1, 0, 1, 1, 1, 1);
1097                break;
1098       case(4): FillReflD0barHistos(d, 0, 0, 0, hypD0bar2, 0, 0, 0, 0);
1099                FillReflD0barHistos(d, 0, 0, 0, hypD0bar2, 1, 1, 1, 1);
1100                break;
1101     }
1102   }
1103   else {
1104     flagLayer1 = 6; flagLayer2 = 6; flagLayer3 = 6; flagLayer4 = 6;
1105     switch(truthHyp) { //fills fhReflD0 and fhReflD0bar with layers containing all true D0 or D0bars, for multi and all hypotheses (layers 2, 0)
1106       case(1): FillReflD0Histos(d, hypD01, 0, 0, 0, 0, 0, 0, 0);
1107                FillReflD0Histos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, 2, 3, 4, 5);
1108                FillReflD0Histos(d, 0, 0, hypD0bar1, hypD0bar2, 2, 3, 6, 6); //merge of opposite particle hyps (D0bar) in layer 6
1109                flagLayer1 = 5;
1110                break;
1111       case(2): FillReflD0Histos(d, 0, hypD02, 0, 0, 0, 0, 0, 0);
1112                FillReflD0Histos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, 3, 2, 4, 5);
1113                FillReflD0Histos(d, 0, 0, hypD0bar1, hypD0bar2, 0, 0, 6, 6); //merge of opposite particle hyps (D0bar) in layer 6
1114                flagLayer2 = 5;
1115                break;
1116       case(3): FillReflD0barHistos(d, 0, 0, hypD0bar1, 0, 0, 0, 0, 0);
1117                FillReflD0barHistos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, 4, 5, 2, 3);
1118                FillReflD0barHistos(d, hypD01, hypD02, 0, 0, 6, 6, 0, 0); //merge of opposite particle hyps (D0) in layer 6
1119                flagLayer3 = 5;
1120                break;
1121       case(4): FillReflD0barHistos(d, 0, 0, 0, hypD0bar2, 0, 0, 0, 0);
1122                FillReflD0barHistos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, 4, 5, 3, 2);
1123                FillReflD0barHistos(d, hypD01, hypD02, 0, 0, 6, 6, 0, 0); //merge of opposite particle hyps (D0) in layer 6
1124                flagLayer4 = 5;
1125                break;
1126     }
1127     FillReflHistos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, flagLayer1, flagLayer2, flagLayer3, flagLayer4); //fill layers 5 and 6 (true and false hyps for multi)
1128     FillReflHistos(d, hypD01, hypD02, hypD0bar1, hypD0bar2, 4, 4, 4, 4); //fill layer 4 (true D0/D0bar - multi hyps only)
1129   }
1130 }
1131
1132 //________________________________________________________________________
1133 Int_t AliAnalysisTaskSESelectHF4Prong::StudyMCTruth(TClonesArray* mcArray, AliAODRecoDecayHF4Prong* d)
1134 {
1135   /*
1136   ---STUDY OF MCTRUTH ON CANDIDATE IN ANALYSIS---
1137   Flag Truth (output):
1138   0 = problems in daughter tracks found
1139   1 = candidate is D01 (piKpipi)
1140   2 = candidate is D02 (pipipiK)
1141   3 = candidate is D0bar1 (Kpipipi)
1142   4 = candidate is D0bar2 (pipiKpi)
1143   */
1144
1145   Int_t truthHyp = 0;
1146
1147   AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
1148   AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1);
1149   AliAODTrack *trk2 = (AliAODTrack*)d->GetDaughter(2);
1150   AliAODTrack *trk3 = (AliAODTrack*)d->GetDaughter(3);
1151   Int_t labels[4];
1152   Int_t pdg[4];
1153   labels[0] = trk0->GetLabel();
1154   labels[1] = trk1->GetLabel();
1155   labels[2] = trk2->GetLabel();
1156   labels[3] = trk3->GetLabel();
1157   if (labels[0]<=0 || labels[1]<=0 || labels[2]<=0 || labels[3]<=0) {AliWarning("Negative Label for daughter, skipping"); return truthHyp;}
1158   AliAODMCParticle* mc0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labels[0]));
1159   AliAODMCParticle* mc1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labels[1]));
1160   AliAODMCParticle* mc2 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labels[2]));
1161   AliAODMCParticle* mc3 = dynamic_cast<AliAODMCParticle*>(mcArray->At(labels[3]));
1162   if (!mc0 || !mc1 || !mc2 || !mc3) {AliWarning("At least one Daughter Particle not found in tree, skipping"); return truthHyp;}
1163   pdg[0] = TMath::Abs(mc0->GetPdgCode());
1164   pdg[1] = TMath::Abs(mc1->GetPdgCode());
1165   pdg[2] = TMath::Abs(mc2->GetPdgCode());
1166   pdg[3] = TMath::Abs(mc3->GetPdgCode());
1167   if (pdg[0]==211 && pdg[1]==321 && pdg[2]==211 && pdg[3]==211) truthHyp = 1;
1168   if (pdg[0]==211 && pdg[1]==211 && pdg[2]==211 && pdg[3]==321) truthHyp = 2;
1169   if (pdg[0]==321 && pdg[1]==211 && pdg[2]==211 && pdg[3]==211) truthHyp = 3;
1170   if (pdg[0]==211 && pdg[1]==211 && pdg[2]==321 && pdg[3]==211) truthHyp = 4;
1171
1172   return truthHyp;
1173 }
1174
1175 //________________________________________________________________________
1176 void AliAnalysisTaskSESelectHF4Prong::FillReflHistos(AliAODRecoDecayHF4Prong* d, Int_t hypD01, Int_t hypD02, Int_t hypD0bar1, Int_t hypD0bar2, Int_t flagLayer1, Int_t flagLayer2, Int_t flagLayer3, Int_t flagLayer4)
1177 {
1178
1179   Double_t ptPart = d->Pt();
1180
1181     if (ptPart > fPtBinH[0]) {
1182       // D01 hyp.
1183       if(hypD01==1) {
1184         if (ptPart < fPtBinH[1])                              fhReflBin1->Fill(fmassD0[0],(double)flagLayer1);
1185         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflBin2->Fill(fmassD0[0],(double)flagLayer1);
1186         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflBin3->Fill(fmassD0[0],(double)flagLayer1);
1187         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflBin4->Fill(fmassD0[0],(double)flagLayer1);
1188         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflBin5->Fill(fmassD0[0],(double)flagLayer1);
1189       }
1190       // D02 hyp.
1191       if(hypD02==1) {
1192         if (ptPart < fPtBinH[1])                              fhReflBin1->Fill(fmassD0[1],(double)flagLayer2);
1193         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflBin2->Fill(fmassD0[1],(double)flagLayer2);
1194         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflBin3->Fill(fmassD0[1],(double)flagLayer2);
1195         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflBin4->Fill(fmassD0[1],(double)flagLayer2);
1196         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflBin5->Fill(fmassD0[1],(double)flagLayer2);
1197       }
1198       // D0bar1 hyp.
1199       if(hypD0bar1==1) {
1200         if (ptPart < fPtBinH[1])                              fhReflBin1->Fill(fmassD0bar[0],(double)flagLayer3);
1201         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflBin2->Fill(fmassD0bar[0],(double)flagLayer3);
1202         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflBin3->Fill(fmassD0bar[0],(double)flagLayer3);
1203         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflBin4->Fill(fmassD0bar[0],(double)flagLayer3);
1204         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflBin5->Fill(fmassD0bar[0],(double)flagLayer3);
1205       }
1206       // D0bar2 hyp.
1207       if(hypD0bar2==1) {
1208         if (ptPart < fPtBinH[1])                              fhReflBin1->Fill(fmassD0bar[1],(double)flagLayer4);
1209         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflBin2->Fill(fmassD0bar[1],(double)flagLayer4);
1210         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflBin3->Fill(fmassD0bar[1],(double)flagLayer4);
1211         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflBin4->Fill(fmassD0bar[1],(double)flagLayer4);
1212         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflBin5->Fill(fmassD0bar[1],(double)flagLayer4);
1213       }
1214     }
1215
1216 }
1217
1218 //________________________________________________________________________
1219 void AliAnalysisTaskSESelectHF4Prong::FillReflD0Histos(AliAODRecoDecayHF4Prong* d, Int_t hypD01, Int_t hypD02, Int_t hypD0bar1, Int_t hypD0bar2, Int_t flagLayforD01, Int_t flagLayforD02, Int_t flagLayforD03, Int_t flagLayforD04)
1220 {
1221
1222   Double_t ptPart = d->Pt();
1223
1224     if (ptPart > fPtBinH[0]) {
1225       // D01 hyp.
1226       if(hypD01==1) {
1227         if (ptPart < fPtBinH[1])                              fhReflD0Bin1->Fill(fmassD0[0],(double)flagLayforD01);
1228         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflD0Bin2->Fill(fmassD0[0],(double)flagLayforD01);
1229         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflD0Bin3->Fill(fmassD0[0],(double)flagLayforD01);
1230         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflD0Bin4->Fill(fmassD0[0],(double)flagLayforD01);
1231         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflD0Bin5->Fill(fmassD0[0],(double)flagLayforD01);
1232       }
1233       // D02 hyp.
1234       if(hypD02==1) {
1235         if (ptPart < fPtBinH[1])                              fhReflD0Bin1->Fill(fmassD0[1],(double)flagLayforD02);
1236         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflD0Bin2->Fill(fmassD0[1],(double)flagLayforD02);
1237         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflD0Bin3->Fill(fmassD0[1],(double)flagLayforD02);
1238         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflD0Bin4->Fill(fmassD0[1],(double)flagLayforD02);
1239         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflD0Bin5->Fill(fmassD0[1],(double)flagLayforD02);
1240       }
1241       // D0bar1 hyp.
1242       if(hypD0bar1==1) {
1243         if (ptPart < fPtBinH[1])                              fhReflD0Bin1->Fill(fmassD0bar[0],(double)flagLayforD03);
1244         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflD0Bin2->Fill(fmassD0bar[0],(double)flagLayforD03);
1245         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflD0Bin3->Fill(fmassD0bar[0],(double)flagLayforD03);
1246         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflD0Bin4->Fill(fmassD0bar[0],(double)flagLayforD03);
1247         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflD0Bin5->Fill(fmassD0bar[0],(double)flagLayforD03);
1248       }
1249       // D0bar2 hyp.
1250       if(hypD0bar2==1) {
1251         if (ptPart < fPtBinH[1])                              fhReflD0Bin1->Fill(fmassD0bar[1],(double)flagLayforD04);
1252         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflD0Bin2->Fill(fmassD0bar[1],(double)flagLayforD04);
1253         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflD0Bin3->Fill(fmassD0bar[1],(double)flagLayforD04);
1254         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflD0Bin4->Fill(fmassD0bar[1],(double)flagLayforD04);
1255         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflD0Bin5->Fill(fmassD0bar[1],(double)flagLayforD04);
1256       }
1257     }
1258
1259 }
1260
1261 //________________________________________________________________________
1262 void AliAnalysisTaskSESelectHF4Prong::FillReflD0barHistos(AliAODRecoDecayHF4Prong* d, Int_t hypD01, Int_t hypD02, Int_t hypD0bar1, Int_t hypD0bar2, Int_t flagLayforD0bar1, Int_t flagLayforD0bar2, Int_t flagLayforD0bar3, Int_t flagLayforD0bar4)
1263 {
1264
1265   Double_t ptPart = d->Pt();
1266
1267     if (ptPart > fPtBinH[0]) {
1268       // D01 hyp.
1269       if(hypD01==1) {
1270         if (ptPart < fPtBinH[1])                              fhReflD0barBin1->Fill(fmassD0[0],(double)flagLayforD0bar1);
1271         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflD0barBin2->Fill(fmassD0[0],(double)flagLayforD0bar1);
1272         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflD0barBin3->Fill(fmassD0[0],(double)flagLayforD0bar1);
1273         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflD0barBin4->Fill(fmassD0[0],(double)flagLayforD0bar1);
1274         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflD0barBin5->Fill(fmassD0[0],(double)flagLayforD0bar1);
1275       }
1276       // D02 hyp.
1277       if(hypD02==1) {
1278         if (ptPart < fPtBinH[1])                              fhReflD0barBin1->Fill(fmassD0[1],(double)flagLayforD0bar2);
1279         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflD0barBin2->Fill(fmassD0[1],(double)flagLayforD0bar2);
1280         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflD0barBin3->Fill(fmassD0[1],(double)flagLayforD0bar2);
1281         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflD0barBin4->Fill(fmassD0[1],(double)flagLayforD0bar2);
1282         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflD0barBin5->Fill(fmassD0[1],(double)flagLayforD0bar2);
1283       }
1284       // D0bar1 hyp.
1285       if(hypD0bar1==1) {
1286         if (ptPart < fPtBinH[1])                              fhReflD0barBin1->Fill(fmassD0bar[0],(double)flagLayforD0bar3);
1287         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflD0barBin2->Fill(fmassD0bar[0],(double)flagLayforD0bar3);
1288         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflD0barBin3->Fill(fmassD0bar[0],(double)flagLayforD0bar3);
1289         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflD0barBin4->Fill(fmassD0bar[0],(double)flagLayforD0bar3);
1290         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflD0barBin5->Fill(fmassD0bar[0],(double)flagLayforD0bar3);
1291       }
1292       // D0bar2 hyp.
1293       if(hypD0bar2==1) {
1294         if (ptPart < fPtBinH[1])                              fhReflD0barBin1->Fill(fmassD0bar[1],(double)flagLayforD0bar4);
1295         else if (ptPart >= fPtBinH[1] && ptPart < fPtBinH[2]) fhReflD0barBin2->Fill(fmassD0bar[1],(double)flagLayforD0bar4);
1296         else if (ptPart >= fPtBinH[2] && ptPart < fPtBinH[3]) fhReflD0barBin3->Fill(fmassD0bar[1],(double)flagLayforD0bar4);
1297         else if (ptPart >= fPtBinH[3] && ptPart < fPtBinH[4]) fhReflD0barBin4->Fill(fmassD0bar[1],(double)flagLayforD0bar4);
1298         else if (ptPart >= fPtBinH[4] && ptPart < fPtBinH[5]) fhReflD0barBin5->Fill(fmassD0bar[1],(double)flagLayforD0bar4);
1299       }
1300     }
1301
1302 }
1303
1304 //________________________________________________________________________
1305 void AliAnalysisTaskSESelectHF4Prong::SetPtBinH(Double_t* ptlimits)
1306 {
1307   for(int i=0; i<6; i++) {fPtBinH[i]=ptlimits[i];}
1308 }
1309
1310 //________________________________________________________________________
1311 void AliAnalysisTaskSESelectHF4Prong::PrintPtBinHandMCFlag()
1312 {
1313   printf("PtBin limits---------\n");
1314   for (int i=0; i<5; i++) {
1315     printf("Bin %d = %.1f to %.1f\n",i+1,fPtBinH[i],fPtBinH[i+1]);
1316   }
1317   printf("MC Truth = %d\n",fMCTruth);
1318   printf("---------------------\n");
1319 }
1320
1321 //________________________________________________________________________
1322 void AliAnalysisTaskSESelectHF4Prong::Terminate(Option_t */*option*/)
1323 {
1324   // Terminate analysis
1325   //
1326
1327 /*
1328   Double_t entries[10] = {0,0,0,0,0,0,0,0,0,0};
1329   for(int i=1;i<=fhReflD0Bin1->GetNbinsX();i++) {
1330     for(int j=1;j<=fhReflD0Bin1->GetNbinsY();j++) {
1331       entries[0] += fhReflD0Bin1->GetBinContent(i,j);
1332       entries[1] += fhReflD0Bin2->GetBinContent(i,j);
1333       entries[2] += fhReflD0Bin3->GetBinContent(i,j);
1334       entries[3] += fhReflD0Bin4->GetBinContent(i,j);
1335       entries[4] += fhReflD0Bin5->GetBinContent(i,j);
1336       entries[5] += fhReflD0barBin1->GetBinContent(i,j);
1337       entries[6] += fhReflD0barBin2->GetBinContent(i,j);
1338       entries[7] += fhReflD0barBin3->GetBinContent(i,j);
1339       entries[8] += fhReflD0barBin4->GetBinContent(i,j);
1340       entries[9] += fhReflD0barBin5->GetBinContent(i,j);
1341     }
1342   }
1343   fhReflD0Bin1->SetEntries(entries[0]);
1344   fhReflD0Bin2->SetEntries(entries[1]);
1345   fhReflD0Bin3->SetEntries(entries[2]);
1346   fhReflD0Bin4->SetEntries(entries[3]);
1347   fhReflD0Bin5->SetEntries(entries[4]);
1348   fhReflD0barBin1->SetEntries(entries[5]);
1349   fhReflD0barBin2->SetEntries(entries[6]);
1350   fhReflD0barBin3->SetEntries(entries[7]);
1351   fhReflD0barBin4->SetEntries(entries[8]);
1352   fhReflD0barBin5->SetEntries(entries[9]);*/
1353
1354   if(fDebug > 1) printf("AnalysisTaskSESelectHF4Prong: Terminate() \n");
1355
1356   fOutput = dynamic_cast<TList*> (GetOutputData(1));
1357   if (!fOutput) {
1358     printf("ERROR: fOutput not available\n");
1359     return;
1360   }
1361   fOutput2 = dynamic_cast<TList*> (GetOutputData(2));
1362   if (!fOutput2) {
1363     printf("ERROR: fOutput not available\n");
1364     return;
1365   }
1366   fOutput3 = dynamic_cast<TList*> (GetOutputData(3));
1367   if (!fOutput3) {
1368     printf("ERROR: fOutput not available\n");
1369     return;
1370   }
1371   fOutput4 = dynamic_cast<TList*> (GetOutputData(4));
1372   if (!fOutput4) {
1373     printf("ERROR: fOutput not available\n");
1374     return;
1375   }
1376   fOutput5 = dynamic_cast<TList*> (GetOutputData(5));
1377   if (!fOutput5) {
1378     printf("ERROR: fOutput not available\n");
1379     return;
1380   }
1381   fOutputC = dynamic_cast<TList*> (GetOutputData(6));
1382   if (!fOutputC) {
1383     printf("ERROR: fOutputC not available\n");
1384     return;
1385   }
1386
1387   fhInvMassD0Sum10MevBin1 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum10MevBin1"));
1388   fhInvMassD0barSum10MevBin1 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum10MevBin1"));
1389   fhInvMassSumAll10MevBin1 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll10MevBin1"));
1390   fhInvMassD0Sum5MevBin1 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum5MevBin1"));
1391   fhInvMassD0barSum5MevBin1 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum5MevBin1"));
1392   fhInvMassSumAll5MevBin1 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll5MevBin1"));
1393   fhInvMassD0Sum10MevBin2 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum10MevBin2"));
1394   fhInvMassD0barSum10MevBin2 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum10MevBin2"));
1395   fhInvMassSumAll10MevBin2 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll10MevBin2"));
1396   fhInvMassD0Sum5MevBin2 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum5MevBin2"));
1397   fhInvMassD0barSum5MevBin2 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum5MevBin2"));
1398   fhInvMassSumAll5MevBin2 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll5MevBin2"));
1399   fhInvMassD0Sum10MevBin3 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum10MevBin3"));
1400   fhInvMassD0barSum10MevBin3 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum10MevBin3"));
1401   fhInvMassSumAll10MevBin3 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll10MevBin3"));
1402   fhInvMassD0Sum5MevBin3 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum5MevBin3"));
1403   fhInvMassD0barSum5MevBin3 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum5MevBin3"));
1404   fhInvMassSumAll5MevBin3 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll5MevBin3"));
1405   fhInvMassD0Sum10MevBin4 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum10MevBin4"));
1406   fhInvMassD0barSum10MevBin4 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum10MevBin4"));
1407   fhInvMassSumAll10MevBin4 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll10MevBin4"));
1408   fhInvMassD0Sum5MevBin4 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum5MevBin4"));
1409   fhInvMassD0barSum5MevBin4 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum5MevBin4"));
1410   fhInvMassSumAll5MevBin4 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll5MevBin4"));
1411   fhInvMassD0Sum10MevBin5 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum10MevBin5"));
1412   fhInvMassD0barSum10MevBin5 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum10MevBin5"));
1413   fhInvMassSumAll10MevBin5 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll10MevBin5"));
1414   fhInvMassD0Sum5MevBin5 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0Sum5MevBin5"));
1415   fhInvMassD0barSum5MevBin5 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassD0barSum5MevBin5"));
1416   fhInvMassSumAll5MevBin5 = dynamic_cast<TH1F*>(fOutput->FindObject("fhInvMassSumAll5MevBin5"));
1417
1418   fScatterP4PID = dynamic_cast<TH2F*>(fOutput->FindObject("fScatterP4PID"));
1419   fPtVsY = dynamic_cast<TH2F*>(fOutputC->FindObject("fPtVsY"));
1420   fPtVsYAll = dynamic_cast<TH2F*>(fOutputC->FindObject("fPtVsYAll"));
1421
1422   fEventCounter = dynamic_cast<TH1F*>(fOutputC->FindObject("fEventCounter"));
1423
1424   fCutDCA = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutDCA"));
1425   fCutDCA3 = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutDCA3"));
1426   fCutDCA2 = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutDCA2"));
1427   fCutDCA5 = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutDCA5"));
1428   fCutVertexDist2 = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutVertexDist2"));
1429   fCutVertexDist3 = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutVertexDist3"));
1430   fCutVertexDist4 = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutVertexDist4"));
1431   fCutCosinePoint = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutCosinePoint"));
1432
1433   fCutPt = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutPt"));
1434   fCutY = dynamic_cast<TH1F*>(fOutputC->FindObject("fCutY"));
1435   fPIDSel = dynamic_cast<TH1F*>(fOutputC->FindObject("fPIDSel"));
1436   fPIDSelBin1 = dynamic_cast<TH1F*>(fOutputC->FindObject("fPIDSelBin1"));
1437   fPIDSelBin2 = dynamic_cast<TH1F*>(fOutputC->FindObject("fPIDSelBin2"));
1438   fPIDSelBin3 = dynamic_cast<TH1F*>(fOutputC->FindObject("fPIDSelBin3"));
1439   fPIDSelBin4 = dynamic_cast<TH1F*>(fOutputC->FindObject("fPIDSelBin4"));
1440   fPIDSelBin5 = dynamic_cast<TH1F*>(fOutputC->FindObject("fPIDSelBin5"));
1441 }
1442