AliPhysics  41af4b0 (41af4b0)
AliAnalysisTaskSEBPlustoD0Pi.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2009, 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  * appeuear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 /* $Id$ */
17 
18 //
19 //
20 // Base class for (BPlus -> D0 pi-> K pi pi) Analysis
21 //
22 //
23 // Cuts are centralized in AliRDHFCutsBPlustoD0Pi
24 // Like sign background is imlemented in the macro
25 //
26 //-----------------------------------------------------------------------
27 //
28 // Author Lennart van Doremalen
29 // Utrecht University - l.v.r.vandoremalen@uu.nl
30 //
31 //-----------------------------------------------------------------------
32 
33 #include <TSystem.h>
34 #include <TChain.h>
35 #include <TParticle.h>
36 #include <TH1I.h>
37 #include "TROOT.h"
38 #include <TDatabasePDG.h>
39 #include <AliAnalysisDataSlot.h>
40 #include <AliAnalysisDataContainer.h>
41 #include "AliRDHFCutsBPlustoD0Pi.h"
42 #include "AliStack.h"
43 #include "AliMCEvent.h"
44 #include "AliAnalysisManager.h"
45 #include "AliAODMCHeader.h"
46 #include "AliAODHandler.h"
47 #include "AliLog.h"
48 #include "AliVertex.h"
49 #include "AliVVertex.h"
50 #include "AliESDVertex.h"
51 #include "AliAODVertex.h"
52 #include "AliVertexerTracks.h"
53 #include "AliExternalTrackParam.h"
54 #include "AliNeutralTrackParam.h"
55 #include "AliAODRecoDecay.h"
56 #include "AliAODRecoDecayHF.h"
59 #include "AliAnalysisVertexingHF.h"
60 #include "AliVertexingHFUtils.h"
61 #include "AliESDtrack.h"
62 #include "AliAODMCParticle.h"
63 #include "AliAODEvent.h"
65 #include "AliAODInputHandler.h"
66 #include <vector>
67 #include <TMatrix.h>
68 #include <TVector3.h>
69 #include <TArrayI.h>
70 #include <bitset>
71 #include <TH3F.h>
72 
73 // #include "TObjectTable.h"
74 
78 
79 //__________________________________________________________________________
82  fListCuts(0),
83  fEvents(0),
84  fUseMCInfo(kFALSE),
85  fOutput(0),
86  fOutputD0FirstDaughter(0),
87  fOutputD0SecondDaughter(0),
88  fOutputBPlusPion(0),
89  fOutputD0(0),
90  fOutputBPlus(0),
91  fOutputD0_D0Pt(0),
92  fOutputBPlusMC(0),
93  fCuts(0),
94  fQuickSignalAnalysis(0),
95  fGetCutInfo(0),
96  fCEvents(0),
97  fBPlusPionTracks(0x0),
98  fD0Tracks(0x0),
99  fShowMask(0),
100  fShowRejection(0),
101  fnPtBins(0),
102  fnPtBinsD0forD0ptbin(0),
103  fnPtBinLimits(0),
104  fnPtBinsD0forD0ptbinLimits(0),
105  fPtBinLimits(0x0),
106  fPtBinLimitsD0forD0ptbin(0x0),
107  fDaughterHistogramArray(),
108  fDaughterHistogramArray2D(),
109  fDaughterHistogramArrayExtra(),
110  fMotherHistogramArray(),
111  fMotherHistogramArray2D(),
112  fMotherHistogramArrayExtra(),
113  fUpgradeSetting(0),
114  fHistMassWindow(0.125),
115  fDegreePerRotation(0),
116  fNumberOfRotations(0),
117  fCheckBackground(0)
118 {
119  //
121  //
122 
123 }
124 //___________________________________________________________________________
126  AliAnalysisTaskSE(name),
127  fListCuts(0),
128  fEvents(0),
129  fUseMCInfo(kFALSE),
130  fOutput(0),
133  fOutputBPlusPion(0),
134  fOutputD0(0),
135  fOutputBPlus(0),
136  fOutputD0_D0Pt(0),
137  fOutputBPlusMC(0),
138  fCuts(0),
140  fGetCutInfo(0),
141  fCEvents(0),
142  fBPlusPionTracks(0x0),
143  fD0Tracks(0x0),
144  fShowMask(0),
145  fShowRejection(0),
146  fnPtBins(0),
148  fnPtBinLimits(0),
150  fPtBinLimits(0x0),
158  fUpgradeSetting(0),
159  fHistMassWindow(0.125),
163 {
164  //
166  //
167 
168  Info("AliAnalysisTaskSEBPlustoD0Pi", "Calling Constructor");
169 
170  fCuts = cuts;
171 
172  DefineInput(0,TChain::Class());
173  DefineOutput(1, TList::Class()); //conters
174  DefineOutput(2, TList::Class()); //My private output
175  DefineOutput(3, TList::Class()); // D0 pion output
176  DefineOutput(4, TList::Class()); // D0 kaon output
177  DefineOutput(5, TList::Class()); // BPlus pion output
178  DefineOutput(6, TList::Class()); // D0 output
179  DefineOutput(7, TList::Class()); // BPlus output
180  DefineOutput(8, TList::Class()); // BPlus output
181  DefineOutput(9, TList::Class()); // BPlusMC output
182 }
183 
184 //___________________________________________________________________________
186  //
188  //
189  Info("~AliAnalysisTaskSEBPlustoD0Pi", "Calling Destructor");
190 
191  delete fOutput;
192  delete fOutputD0FirstDaughter;
194  delete fOutputBPlusPion;
195  delete fOutputD0;
196  delete fOutputBPlus;
197  delete fOutputD0_D0Pt;
198  delete fOutputBPlusMC;
199  delete fCuts;
200  delete fCEvents;
201  delete fD0Tracks;
202  delete fBPlusPionTracks;
203  delete fListCuts;
204  delete fPtBinLimits;
206 
207  for (Int_t i=0;i<4;i++) {
208  for (Int_t j=0;j<6;j++) {
209  for (Int_t k=0;k<15;k++) {
210  delete fDaughterHistogramArray[i][j][k]; fDaughterHistogramArray[i][j][k] = nullptr;
211  }
212  }
213  }
214  for (Int_t i=0;i<4;i++) {
215  for (Int_t j=0;j<6;j++) {
216  delete fDaughterHistogramArray2D[i][j]; fDaughterHistogramArray2D[i][j] = nullptr;
217  }
218  }
219  for (Int_t i=0;i<4;i++) {
220  for (Int_t j=0;j<6;j++) {
221  delete fDaughterHistogramArrayExtra[i][j]; fDaughterHistogramArrayExtra[i][j] = nullptr;
222  }
223  }
224  for (Int_t i=0;i<6;i++) {
225  for (Int_t j=0;j<99;j++) {
226  for (Int_t k=0;k<60;k++) {
227  delete fMotherHistogramArray[i][j][k]; fMotherHistogramArray[i][j][k] = nullptr;
228  }
229  }
230  }
231  for (Int_t i=0;i<6;i++) {
232  for (Int_t j=0;j<99;j++) {
233  for (Int_t k=0;k<60;k++) {
234  delete fMotherHistogramArray2D[i][j][k]; fMotherHistogramArray2D[i][j][k] = nullptr;
235  }
236  }
237  }
238  for (Int_t i=0;i<7;i++) {
239  for (Int_t j=0;j<10;j++) {
240  delete fMotherHistogramArrayExtra[i][j]; fMotherHistogramArrayExtra[i][j] = nullptr;
241  }
242  }
243 }
244 //_________________________________________________
246  //
248  //
249 
250  if (fDebug > 1) printf("AliAnalysisTaskSEBPlustoD0Pi::Init() \n");
251 
252  return;
253 }
254 //_________________________________________________
256 
257  //==================================================================================
258  // USER EXECUTION FUNCTION - start
259  //==================================================================================
260  //
261  // This is the main function that has to be altered to do heavy flavour analysis.
262  //
263  //==================================================================================
264 
265  if (!fInputEvent) {
266  Error("UserExec", "NO EVENT FOUND!");
267  return;
268  }
269 
270  if (fEvents % 50 == 0) {
271  std::cout << "\r" << "Analysing event number: " << fEvents << std::endl;
272  }
273 
274  fEvents++;
275 
276  // Show trigger mask
277  if (fShowMask)
278  {
279  std::bitset<32> maskEV(((AliAODInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
280  std::cout << "Event mask: " << maskEV << std::endl;
281  std::cout << "Trigger mask: " << std::bitset<32>(fCuts->GetTriggerMask()) << std::endl;
282  }
283 
284 
285  //==================================================================================
286  // EVENT INITIALIZATION - start
287  //==================================================================================
288 
289 
290  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
291  TClonesArray * D0TracksFromFriendFile = 0;
292  fCEvents->Fill(1);
293 
294  if (!aodEvent && AODEvent() && IsStandardAOD()) {
295  // In case there is an AOD handler writing a standard AOD, use the AOD
296  // event in memory rather than the input (ESD) event.
297  aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
298  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
299  // have to taken from the AOD event hold by the AliAODExtension
300  AliAODHandler* aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
301  if (aodHandler->GetExtensions()) {
302  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
303  AliAODEvent *aodFromExt = ext->GetAOD();
304  D0TracksFromFriendFile = (TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
305  }
306  } else {
307  D0TracksFromFriendFile = (TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
308  }
309 
310  // fix for temporary bug in ESDfilter
311  // the AODs with null vertex pointer didn't pass the PhysSel
312  if (!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField()) < 0.001) return;
313  fCEvents->Fill(2);
314 
315  // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
316  TString trigclass = aodEvent->GetFiredTriggerClasses();
317  if (trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fCEvents->Fill(5);
318 
319  if(!fCuts->IsEventSelected(aodEvent))
320  {
321  if(fShowRejection) std::cout << "Event rejected by code: " << fCuts->GetWhyRejection() << std::endl;
322  if(fCuts->GetWhyRejection()==6) // rejected for Z vertex
323  {
324  fCEvents->Fill(6);
325  return;
326  }
327  }
328 
329  Bool_t isEvSel = fCuts->IsEventSelected(aodEvent);
330  if (!isEvSel) return;
331  fCEvents->Fill(3);
332 
333  //get the magnetic field
334  Double_t bz = (Double_t)aodEvent->GetMagneticField();
335 
336  // AOD primary vertex
337  AliAODVertex *primaryVertex = (AliAODVertex*)aodEvent->GetPrimaryVertex();
338  if (!primaryVertex) return;
339  if (primaryVertex->GetNContributors() < 1) return;
340  fCEvents->Fill(4);
341 
342  if(!D0TracksFromFriendFile)
343  {
344  AliInfo("Could not find array of HF vertices, skipping the event");
345  return;
346  }
347  else AliDebug(2, Form("Found %d vertices",D0TracksFromFriendFile->GetEntriesFast()));
348 
349  AliAODMCHeader *mcHeader = 0;
350  if(fUseMCInfo)
351  {
352  mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
353  if(!mcHeader) {
354  printf(" MC header branch not found!\n");
355  return;
356  }
357  }
358 
359 
360  //==================================================================================
361  // EVENT INITIALIZATION - end
362  //==================================================================================
363  // BPlus LOSS BY CUTS TRACKER - start
364  //==================================================================================
365 
366  TClonesArray *mcTrackArray = nullptr;
367  if(fUseMCInfo) mcTrackArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
368  if(fUseMCInfo && !mcTrackArray) return;
369 
370  //Here we create an array and vectors that we use to look how many BPluss can be reconstructed after each cut
371  TMatrix * BPlustoD0PiLabelMatrix = new TMatrix(0, 5);
372 
373  //we fill the array with all BPlus->D0Pi tracks
374  if (fUseMCInfo) {
375  BPlustoD0PiSignalTracksInMC(mcTrackArray, aodEvent, BPlustoD0PiLabelMatrix, fOutputBPlusMC);
376  }
377 
378  //==================================================================================
379  // BPlus LOSS BY CUTS TRACKER - end
380  //==================================================================================
381  // PARTICLE SELECTION LOOP - start
382  //==================================================================================
383  //
384  // Here we select and reconstruct the particles for the BPlus->D*Pion decay.
385  //
386  //==================================================================================
387 
388 
389  BPlusPionSelection(aodEvent, primaryVertex, bz, mcTrackArray, BPlustoD0PiLabelMatrix, mcHeader);
390  D0Selection(aodEvent, primaryVertex, bz, mcTrackArray, BPlustoD0PiLabelMatrix, D0TracksFromFriendFile, mcHeader);
391  BPlusSelection(aodEvent, primaryVertex, bz, mcTrackArray, BPlustoD0PiLabelMatrix, D0TracksFromFriendFile, mcHeader);
392 
393  // Clear arrays and memory management:
394  fD0Tracks->erase(fD0Tracks->begin(), fD0Tracks->end());
395  fBPlusPionTracks->erase(fBPlusPionTracks->begin(), fBPlusPionTracks->end());
396 
397  delete BPlustoD0PiLabelMatrix; BPlustoD0PiLabelMatrix = NULL;
398 
399  //==================================================================================
400  // PARTICLE SELECTION LOOP - end
401  //==================================================================================
402 
403  PostData(1, fOutput);
404  PostData(3, fOutputD0FirstDaughter);
405  PostData(4, fOutputD0SecondDaughter);
406  PostData(5, fOutputBPlusPion);
407  PostData(6, fOutputD0);
408  PostData(7, fOutputBPlus);
409  PostData(8, fOutputD0_D0Pt);
410  PostData(9, fOutputBPlusMC);
411 
412 
413  //==================================================================================
414  // USER EXECUTION FUNCTION - end
415  //==================================================================================
416 
417 }
418 //___________________________________ terminate ___________________________
423 
424  //Info("Terminate","");
425  AliAnalysisTaskSE::Terminate();
426 
427  fOutput = dynamic_cast<TList*> (GetOutputData(1));
428  if (!fOutput) {
429  printf("ERROR: fOutput not available\n");
430  return;
431  }
432 
433  fCEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fCEvents"));
434 
435  fListCuts = dynamic_cast<TList*> (GetOutputData(2));
436  if (!fListCuts) {
437  printf("ERROR: fListCuts not available\n");
438  return;
439  }
440  fOutputD0FirstDaughter = dynamic_cast<TList*> (GetOutputData(3));
441  if (!fOutputD0FirstDaughter) {
442  printf("ERROR: fOutputD0FirstDaughter not available\n");
443  return;
444  }
445  fOutputD0SecondDaughter = dynamic_cast<TList*> (GetOutputData(4));
447  printf("ERROR: fOutputD0SecondDaughter not available\n");
448  return;
449  }
450  fOutputBPlusPion = dynamic_cast<TList*> (GetOutputData(5));
451  if (!fOutputBPlusPion) {
452  printf("ERROR: fOutputBPlusPion not available\n");
453  return;
454  }
455  fOutputD0 = dynamic_cast<TList*> (GetOutputData(6));
456  if (!fOutputD0) {
457  printf("ERROR: fOutputD0 not available\n");
458  return;
459  }
460  fOutputBPlus = dynamic_cast<TList*> (GetOutputData(7));
461  if (!fOutputBPlus) {
462  printf("ERROR: fOutputBPlus not available\n");
463  return;
464  }
465  fOutputD0_D0Pt = dynamic_cast<TList*> (GetOutputData(8));
466  if (!fOutputD0_D0Pt) {
467  printf("ERROR: fOutputD0_D0Pt not available\n");
468  return;
469  }
470  fOutputBPlusMC = dynamic_cast<TList*> (GetOutputData(9));
471  if (!fOutputBPlusMC) {
472  printf("ERROR: fOutputBPlusMC not available\n");
473  return;
474  }
475  return;
476 }
477 
478 //___________________________________________________________________________
481  Info("UserCreateOutputObjects", "CreateOutputObjects of task %s\n", GetName());
482 
483  //slot #1
484  //OpenFile(1);
485  fOutput = new TList();
486  fOutput->SetOwner();
487  fOutput->SetName("chist0");
488 
490  fOutputD0FirstDaughter->SetOwner();
491  fOutputD0FirstDaughter->SetName("listD0Pion");
492 
494  fOutputD0SecondDaughter->SetOwner();
495  fOutputD0SecondDaughter->SetName("listD0Kaon");
496 
497  fOutputBPlusPion = new TList();
498  fOutputBPlusPion->SetOwner();
499  fOutputBPlusPion->SetName("listBPlusPion");
500 
501  fOutputD0 = new TList();
502  fOutputD0->SetOwner();
503  fOutputD0->SetName("listD0");
504 
505  fOutputBPlus = new TList();
506  fOutputBPlus->SetOwner();
507  fOutputBPlus->SetName("listBPlus");
508 
509  fOutputD0_D0Pt = new TList();
510  fOutputD0_D0Pt->SetOwner();
511  fOutputD0_D0Pt->SetName("listD0_D0Pt");
512 
513  fOutputBPlusMC = new TList();
514  fOutputBPlusMC->SetOwner();
515  fOutputBPlusMC->SetName("listBPlusMC");
516 
517  // we prepare vectors that will save the positions of the daughter tracks in the track list during the reconstruction
518  fBPlusPionTracks = new std::vector<Int_t>;
519  fD0Tracks = new std::vector<Int_t>;
520 
521  // we get information on the pt bins
523  fnPtBinLimits = fnPtBins + 1;
525 
529 
530 
531  std::cout << "Nr. of BPlus meson bins: " << fCuts->GetNPtBins() << " limits: " << std::endl;
532  for (int i = 0; i < fnPtBinLimits; ++i)
533  {
534  std::cout << fPtBinLimits[i] << " " << std::endl;
535  }
536  std::cout << std::endl;
537  std::cout << "Nr. of D0 meson bins: " << fCuts->GetNPtBinsD0forD0ptbin() << " limits: " << std::endl;
538  for (int i = 0; i < fnPtBinsD0forD0ptbinLimits; ++i)
539  {
540  std::cout << fPtBinLimitsD0forD0ptbin[i] << " " << std::endl;
541  }
542  std::cout << std::endl;
543 
544 
545  fListCuts=new TList();
546  fListCuts->SetOwner();
547  fListCuts->SetName("Cuts");
549  // Post the data
550  fListCuts->Add(copyfCuts);
551 
552  // define histograms
554 
555  PostData(1, fOutput);
556  PostData(2, fListCuts);
557  PostData(3, fOutputD0FirstDaughter);
558  PostData(4, fOutputD0SecondDaughter);
559  PostData(5, fOutputBPlusPion);
560  PostData(6, fOutputD0);
561  PostData(7, fOutputBPlus);
562  PostData(8, fOutputD0_D0Pt);
563  PostData(9, fOutputBPlusMC);
564 
565  return;
566 }
567 //___________________________________ histograms _______________________________________
570 
571 
572  fCEvents = new TH1F("fCEvents","conter",13,0,13);
573  fCEvents->SetStats(kTRUE);
574  fCEvents->GetXaxis()->SetTitle("1");
575  fCEvents->GetYaxis()->SetTitle("counts");
576  fCEvents->GetXaxis()->SetBinLabel(2,"no. of events");
577  fCEvents->GetXaxis()->SetBinLabel(3,"good prim vtx and B field");
578  fCEvents->GetXaxis()->SetBinLabel(4,"no event selected");
579  fCEvents->GetXaxis()->SetBinLabel(5,"no vtx contributors");
580  fCEvents->GetXaxis()->SetBinLabel(6,"trigger for PbPb");
581  fCEvents->GetXaxis()->SetBinLabel(7,"no z vtx");
582  fCEvents->GetXaxis()->SetBinLabel(8,"...");
583  fCEvents->GetXaxis()->SetBinLabel(12,"no. of D0 fail to be rec");
584  fOutput->Add(fCEvents);
585 
586  //====================================================
587 
588  TString name_mc_BPlus_pt ="mc_BPlus_pt";
589  TH1F* hist_mc_BPlus_pt = new TH1F(name_mc_BPlus_pt.Data(),"Pt monte carlo BPlus in BPlus->D*#pi; p_{T} [GeV/c]; Entries",400,0,20);
590  hist_mc_BPlus_pt->Sumw2();
591  hist_mc_BPlus_pt->SetLineColor(6);
592  hist_mc_BPlus_pt->SetMarkerStyle(20);
593  hist_mc_BPlus_pt->SetMarkerSize(0.6);
594  hist_mc_BPlus_pt->SetMarkerColor(6);
595  TH1F* histogram_mc_BPlus_pt = (TH1F*)hist_mc_BPlus_pt->Clone();
596  fOutputBPlusMC->Add(histogram_mc_BPlus_pt);
597 
598  TString name_mc_BPlus_pion_pt ="mc_BPlus_pion_pt";
599  TH1F* hist_mc_BPlus_pion_pt = new TH1F(name_mc_BPlus_pion_pt.Data(),"Pt monte carlo pion of BPlus in BPlus->D*#pi; p_{T} [GeV/c]; Entries",400,0,20);
600  hist_mc_BPlus_pion_pt->Sumw2();
601  hist_mc_BPlus_pion_pt->SetLineColor(6);
602  hist_mc_BPlus_pion_pt->SetMarkerStyle(20);
603  hist_mc_BPlus_pion_pt->SetMarkerSize(0.6);
604  hist_mc_BPlus_pion_pt->SetMarkerColor(6);
605  TH1F* histogram_mc_BPlus_pion_pt = (TH1F*)hist_mc_BPlus_pion_pt->Clone();
606  fOutputBPlusMC->Add(histogram_mc_BPlus_pion_pt);
607 
608  TString name_mc_D0_pt ="mc_D0_pt";
609  TH1F* hist_mc_D0_pt = new TH1F(name_mc_D0_pt.Data(),"Pt monte carlo D0 in BPlus->D*#pi; p_{T} [GeV/c]; Entries",400,0,20);
610  hist_mc_D0_pt->Sumw2();
611  hist_mc_D0_pt->SetLineColor(6);
612  hist_mc_D0_pt->SetMarkerStyle(20);
613  hist_mc_D0_pt->SetMarkerSize(0.6);
614  hist_mc_D0_pt->SetMarkerColor(6);
615  TH1F* histogram_mc_D0_pt = (TH1F*)hist_mc_D0_pt->Clone();
616  fOutputBPlusMC->Add(histogram_mc_D0_pt);
617 
618  TString name_mc_D0_pion_pt ="mc_D0_pion_pt";
619  TH1F* hist_mc_D0_pion_pt = new TH1F(name_mc_D0_pion_pt.Data(),"Pt monte carlo pion of D0 in BPlus->D*#pi; p_{T} [GeV/c]; Entries",400,0,20);
620  hist_mc_D0_pion_pt->Sumw2();
621  hist_mc_D0_pion_pt->SetLineColor(6);
622  hist_mc_D0_pion_pt->SetMarkerStyle(20);
623  hist_mc_D0_pion_pt->SetMarkerSize(0.6);
624  hist_mc_D0_pion_pt->SetMarkerColor(6);
625  TH1F* histogram_mc_D0_pion_pt = (TH1F*)hist_mc_D0_pion_pt->Clone();
626  fOutputBPlusMC->Add(histogram_mc_D0_pion_pt);
627 
628  TString name_mc_D0_kaon_pt ="mc_D0_kaon_pt";
629  TH1F* hist_mc_D0_kaon_pt = new TH1F(name_mc_D0_kaon_pt.Data(),"Pt monte carlo kaon of D0 in BPlus->D*#pi; p_{T} [GeV/c]; Entries",400,0,20);
630  hist_mc_D0_kaon_pt->Sumw2();
631  hist_mc_D0_kaon_pt->SetLineColor(6);
632  hist_mc_D0_kaon_pt->SetMarkerStyle(20);
633  hist_mc_D0_kaon_pt->SetMarkerSize(0.6);
634  hist_mc_D0_kaon_pt->SetMarkerColor(6);
635  TH1F* histogram_mc_D0_kaon_pt = (TH1F*)hist_mc_D0_kaon_pt->Clone();
636  fOutputBPlusMC->Add(histogram_mc_D0_kaon_pt);
637 
638  TString name_mc_BPlus_rapidity_true ="mc_BPlus_rapidity_true";
639  TH1F* hist_mc_BPlus_rapidity_true = new TH1F(name_mc_BPlus_rapidity_true.Data(),"rapidity_true monte carlo BPlus in BPlus->D*#pi; Y; Entries",5000,-20,20);
640  hist_mc_BPlus_rapidity_true->Sumw2();
641  hist_mc_BPlus_rapidity_true->SetLineColor(6);
642  hist_mc_BPlus_rapidity_true->SetMarkerStyle(20);
643  hist_mc_BPlus_rapidity_true->SetMarkerSize(0.6);
644  hist_mc_BPlus_rapidity_true->SetMarkerColor(6);
645  TH1F* histogram_mc_BPlus_rapidity_true = (TH1F*)hist_mc_BPlus_rapidity_true->Clone();
646  fOutputBPlusMC->Add(histogram_mc_BPlus_rapidity_true);
647 
648  TString name_mc_BPlus_pion_rapidity_true ="mc_BPlus_pion_rapidity_true";
649  TH1F* hist_mc_BPlus_pion_rapidity_true = new TH1F(name_mc_BPlus_pion_rapidity_true.Data(),"rapidity_true monte carlo pion of BPlus in BPlus->D*#pi; Y; Entries",5000,-20,20);
650  hist_mc_BPlus_pion_rapidity_true->Sumw2();
651  hist_mc_BPlus_pion_rapidity_true->SetLineColor(6);
652  hist_mc_BPlus_pion_rapidity_true->SetMarkerStyle(20);
653  hist_mc_BPlus_pion_rapidity_true->SetMarkerSize(0.6);
654  hist_mc_BPlus_pion_rapidity_true->SetMarkerColor(6);
655  TH1F* histogram_mc_BPlus_pion_rapidity_true = (TH1F*)hist_mc_BPlus_pion_rapidity_true->Clone();
656  fOutputBPlusMC->Add(histogram_mc_BPlus_pion_rapidity_true);
657 
658  TString name_mc_D0_rapidity_true ="mc_D0_rapidity_true";
659  TH1F* hist_mc_D0_rapidity_true = new TH1F(name_mc_D0_rapidity_true.Data(),"rapidity_true monte carlo D0 in BPlus->D*#pi; Y; Entries",5000,-20,20);
660  hist_mc_D0_rapidity_true->Sumw2();
661  hist_mc_D0_rapidity_true->SetLineColor(6);
662  hist_mc_D0_rapidity_true->SetMarkerStyle(20);
663  hist_mc_D0_rapidity_true->SetMarkerSize(0.6);
664  hist_mc_D0_rapidity_true->SetMarkerColor(6);
665  TH1F* histogram_mc_D0_rapidity_true = (TH1F*)hist_mc_D0_rapidity_true->Clone();
666  fOutputBPlusMC->Add(histogram_mc_D0_rapidity_true);
667 
668  TString name_mc_D0_pion_rapidity_true ="mc_D0_pion_rapidity_true";
669  TH1F* hist_mc_D0_pion_rapidity_true = new TH1F(name_mc_D0_pion_rapidity_true.Data(),"rapidity_true monte carlo pion of D0 in BPlus->D*#pi; Y; Entries",5000,-20,20);
670  hist_mc_D0_pion_rapidity_true->Sumw2();
671  hist_mc_D0_pion_rapidity_true->SetLineColor(6);
672  hist_mc_D0_pion_rapidity_true->SetMarkerStyle(20);
673  hist_mc_D0_pion_rapidity_true->SetMarkerSize(0.6);
674  hist_mc_D0_pion_rapidity_true->SetMarkerColor(6);
675  TH1F* histogram_mc_D0_pion_rapidity_true = (TH1F*)hist_mc_D0_pion_rapidity_true->Clone();
676  fOutputBPlusMC->Add(histogram_mc_D0_pion_rapidity_true);
677 
678  TString name_mc_D0_kaon_rapidity_true ="mc_D0_kaon_rapidity_true";
679  TH1F* hist_mc_D0_kaon_rapidity_true = new TH1F(name_mc_D0_kaon_rapidity_true.Data(),"rapidity_true monte carlo kaon of D0 in BPlus->D*#pi; Y; Entries",5000,-20,20);
680  hist_mc_D0_kaon_rapidity_true->Sumw2();
681  hist_mc_D0_kaon_rapidity_true->SetLineColor(6);
682  hist_mc_D0_kaon_rapidity_true->SetMarkerStyle(20);
683  hist_mc_D0_kaon_rapidity_true->SetMarkerSize(0.6);
684  hist_mc_D0_kaon_rapidity_true->SetMarkerColor(6);
685  TH1F* histogram_mc_D0_kaon_rapidity_true = (TH1F*)hist_mc_D0_kaon_rapidity_true->Clone();
686  fOutputBPlusMC->Add(histogram_mc_D0_kaon_rapidity_true);
687 
688  TString name_mc_BPlus_pseudorapidity_true ="mc_BPlus_pseudorapidity_true";
689  TH1F* hist_mc_BPlus_pseudorapidity_true = new TH1F(name_mc_BPlus_pseudorapidity_true.Data(),"pseudorapidity_true monte carlo BPlus in BPlus->D*#pi; #eta; Entries",5000,-20,20);
690  hist_mc_BPlus_pseudorapidity_true->Sumw2();
691  hist_mc_BPlus_pseudorapidity_true->SetLineColor(6);
692  hist_mc_BPlus_pseudorapidity_true->SetMarkerStyle(20);
693  hist_mc_BPlus_pseudorapidity_true->SetMarkerSize(0.6);
694  hist_mc_BPlus_pseudorapidity_true->SetMarkerColor(6);
695  TH1F* histogram_mc_BPlus_pseudorapidity_true = (TH1F*)hist_mc_BPlus_pseudorapidity_true->Clone();
696  fOutputBPlusMC->Add(histogram_mc_BPlus_pseudorapidity_true);
697 
698  TString name_mc_BPlus_pion_pseudorapidity_true ="mc_BPlus_pion_pseudorapidity_true";
699  TH1F* hist_mc_BPlus_pion_pseudorapidity_true = new TH1F(name_mc_BPlus_pion_pseudorapidity_true.Data(),"pseudorapidity_true monte carlo pion of BPlus in BPlus->D*#pi; #eta; Entries",5000,-20,20);
700  hist_mc_BPlus_pion_pseudorapidity_true->Sumw2();
701  hist_mc_BPlus_pion_pseudorapidity_true->SetLineColor(6);
702  hist_mc_BPlus_pion_pseudorapidity_true->SetMarkerStyle(20);
703  hist_mc_BPlus_pion_pseudorapidity_true->SetMarkerSize(0.6);
704  hist_mc_BPlus_pion_pseudorapidity_true->SetMarkerColor(6);
705  TH1F* histogram_mc_BPlus_pion_pseudorapidity_true = (TH1F*)hist_mc_BPlus_pion_pseudorapidity_true->Clone();
706  fOutputBPlusMC->Add(histogram_mc_BPlus_pion_pseudorapidity_true);
707 
708  TString name_mc_D0_pseudorapidity_true ="mc_D0_pseudorapidity_true";
709  TH1F* hist_mc_D0_pseudorapidity_true = new TH1F(name_mc_D0_pseudorapidity_true.Data(),"pseudorapidity_true monte carlo D0 in BPlus->D*#pi; #eta; Entries",5000,-20,20);
710  hist_mc_D0_pseudorapidity_true->Sumw2();
711  hist_mc_D0_pseudorapidity_true->SetLineColor(6);
712  hist_mc_D0_pseudorapidity_true->SetMarkerStyle(20);
713  hist_mc_D0_pseudorapidity_true->SetMarkerSize(0.6);
714  hist_mc_D0_pseudorapidity_true->SetMarkerColor(6);
715  TH1F* histogram_mc_D0_pseudorapidity_true = (TH1F*)hist_mc_D0_pseudorapidity_true->Clone();
716  fOutputBPlusMC->Add(histogram_mc_D0_pseudorapidity_true);
717 
718  TString name_mc_D0_pion_pseudorapidity_true ="mc_D0_pion_pseudorapidity_true";
719  TH1F* hist_mc_D0_pion_pseudorapidity_true = new TH1F(name_mc_D0_pion_pseudorapidity_true.Data(),"pseudorapidity_true monte carlo pion of D0 in BPlus->D*#pi; #eta; Entries",5000,-20,20);
720  hist_mc_D0_pion_pseudorapidity_true->Sumw2();
721  hist_mc_D0_pion_pseudorapidity_true->SetLineColor(6);
722  hist_mc_D0_pion_pseudorapidity_true->SetMarkerStyle(20);
723  hist_mc_D0_pion_pseudorapidity_true->SetMarkerSize(0.6);
724  hist_mc_D0_pion_pseudorapidity_true->SetMarkerColor(6);
725  TH1F* histogram_mc_D0_pion_pseudorapidity_true = (TH1F*)hist_mc_D0_pion_pseudorapidity_true->Clone();
726  fOutputBPlusMC->Add(histogram_mc_D0_pion_pseudorapidity_true);
727 
728  TString name_mc_D0_kaon_pseudorapidity_true ="mc_D0_kaon_pseudorapidity_true";
729  TH1F* hist_mc_D0_kaon_pseudorapidity_true = new TH1F(name_mc_D0_kaon_pseudorapidity_true.Data(),"pseudorapidity_true monte carlo kaon of D0 in BPlus->D*#pi; #eta; Entries",5000,-20,20);
730  hist_mc_D0_kaon_pseudorapidity_true->Sumw2();
731  hist_mc_D0_kaon_pseudorapidity_true->SetLineColor(6);
732  hist_mc_D0_kaon_pseudorapidity_true->SetMarkerStyle(20);
733  hist_mc_D0_kaon_pseudorapidity_true->SetMarkerSize(0.6);
734  hist_mc_D0_kaon_pseudorapidity_true->SetMarkerColor(6);
735  TH1F* histogram_mc_D0_kaon_pseudorapidity_true = (TH1F*)hist_mc_D0_kaon_pseudorapidity_true->Clone();
736  fOutputBPlusMC->Add(histogram_mc_D0_kaon_pseudorapidity_true);
737 
738 
739  //==================================================
740 
741  TString name_BPluss_in_analysis ="BPlus_in_analysis";
742  TH1F* hist_BPluss_in_analysis = new TH1F(name_BPluss_in_analysis.Data(),"Number of BPluss to kpipipi in the Analysis; Entries",10,0,10);
743  hist_BPluss_in_analysis->Sumw2();
744  hist_BPluss_in_analysis->SetLineColor(6);
745  hist_BPluss_in_analysis->SetMarkerStyle(20);
746  hist_BPluss_in_analysis->SetMarkerSize(0.6);
747  hist_BPluss_in_analysis->SetMarkerColor(6);
748  hist_BPluss_in_analysis->SetStats(kTRUE);
749  hist_BPluss_in_analysis->GetXaxis()->SetBinLabel(1,"no. of BPlus");
750  hist_BPluss_in_analysis->GetXaxis()->SetBinLabel(2,"no. of BPlus to kpipi");
751  hist_BPluss_in_analysis->GetXaxis()->SetBinLabel(3,"no. with all tracks in event");
752  hist_BPluss_in_analysis->GetXaxis()->SetBinLabel(4,"no. ...");
753  hist_BPluss_in_analysis->GetXaxis()->SetBinLabel(5,"no. ...");
754  hist_BPluss_in_analysis->GetXaxis()->SetBinLabel(6,"no. ...");
755  hist_BPluss_in_analysis->GetXaxis()->SetBinLabel(7,"no. ...");
756  hist_BPluss_in_analysis->GetXaxis()->SetBinLabel(8,"no. ...");
757  TH1F* hist_BPluss_in_analysis_mc = (TH1F*)hist_BPluss_in_analysis->Clone();
758  fOutputBPlusMC->Add(hist_BPluss_in_analysis_mc);
759 
760  TString name_BPlus_per_bin ="BPlus_per_bin";
761  TH1F* hist_BPlus_per_bin = new TH1F(name_BPlus_per_bin.Data(),"Number of BPlus to kpipi in the Analysis per bin; Entries",fnPtBins,0,fnPtBins);
762  for (Int_t i = 0; i < fnPtBins; ++i)
763  {
764  TString bin_name = "";
765  bin_name += fPtBinLimits[i];
766  bin_name += "-";
767  bin_name += fPtBinLimits[i+1];
768  hist_BPlus_per_bin->GetXaxis()->SetBinLabel(i+1,bin_name);
769  }
770  TH1F* hist_BPlus_per_bin_mc = (TH1F*)hist_BPlus_per_bin->Clone();
771  fOutputBPlusMC->Add(hist_BPlus_per_bin_mc);
772 
773  TString name_BPlus_per_bin_in_Acc ="BPlus_per_bin_in_Acc";
774  TH1F* hist_BPlus_per_bin_in_Acc = new TH1F(name_BPlus_per_bin_in_Acc.Data(),"Number of BPlus to kpipi in the Analysis per bin with all daughters in acceptance; Entries",fnPtBins,0,fnPtBins);
775  for (Int_t i = 0; i < fnPtBins; ++i)
776  {
777  TString bin_name = "";
778  bin_name += fPtBinLimits[i];
779  bin_name += "-";
780  bin_name += fPtBinLimits[i+1];
781  hist_BPlus_per_bin_in_Acc->GetXaxis()->SetBinLabel(i+1,bin_name);
782  }
783  TH1F* hist_BPlus_per_bin_in_Acc_mc = (TH1F*)hist_BPlus_per_bin_in_Acc->Clone();
784  fOutputBPlusMC->Add(hist_BPlus_per_bin_in_Acc_mc);
785 
786  //======================================================================================================================================================
787 
788  //we make the histograms for the Pions and Kaon
789  for (Int_t i = 0; i < 3; i++){
790 
791  TString add_name = "";
792  TList * listout;
793  if(i==0) listout = fOutputD0FirstDaughter;
794  if(i==1) listout = fOutputD0SecondDaughter;
795  if(i==2) listout = fOutputBPlusPion;
796 
797  for (Int_t j = 0; j < 6; j++){
798  if(j==0) add_name = "";
799  if(j==1) add_name = "Signal";
800  if(j==2) add_name = "Cut";
801  if(j==3) add_name = "SignalCut";
802  if(j==4) add_name = "Result";
803  if(j==5) add_name = "SignalResult";
804 
805  TString name_Histogram = "";
806  TString discription_Histogram = "";
807  Int_t numberOfBins = 0;
808  Double_t lowerBound = 0.0;
809  Double_t upperBound = 0.0;
810 
811  for (Int_t k = 0; k < 10; ++k)
812  {
813  if(k==0){name_Histogram = "ptTrack"; discription_Histogram = "pt track; p_{T} [GeV/c]; Entries"; numberOfBins = 600; lowerBound = 0; upperBound = 30;}
814  if(k==1){name_Histogram = "momentumTrack"; discription_Histogram = "momentum track; p [GeV/c]; Entries"; numberOfBins = 600; lowerBound = 0; upperBound = 30;}
815  if(k==2){name_Histogram = "numberOfITS"; discription_Histogram = "Number of ITS clusters track; [#]; Entries"; numberOfBins = 10; lowerBound = -0.5; upperBound = 9.5;}
816  if(k==3){name_Histogram = "numberOfTPC"; discription_Histogram = "Number of TPC clusters track; [#]; Entries"; numberOfBins = 601; lowerBound = -0.5; upperBound = 600.5;}
817  if(k==4){name_Histogram = "pointsOnITS"; discription_Histogram = "Number of ITS clusters track per layer; [#]; Entries"; numberOfBins = 10; lowerBound = -0.5; upperBound = 9.5;}
818  if(k==5){name_Histogram = "nSigmaTPC"; discription_Histogram = "n sigma TPC for track PID; sigma; Entries"; numberOfBins = 500; lowerBound = -5; upperBound = 5;}
819  if(k==6){name_Histogram = "nSigmaTOF"; discription_Histogram = "n sigma TOF for track PID; sigma; Entries"; numberOfBins = 500; lowerBound = -5; upperBound = 5;}
820  if(k==7){name_Histogram = "nSigmaTPCandTOF"; discription_Histogram = "n sigma TPC and TOF for track PID; a.u.; Entries"; numberOfBins = 1000; lowerBound = 0; upperBound = 10;}
821  if(k==8){name_Histogram = "impactParameter"; discription_Histogram = "Impact Parameter track; [cm]; Entries"; numberOfBins = 2000; lowerBound = 0; upperBound = 0.5;}
822  if(k==9){name_Histogram = "EtaTrack"; discription_Histogram = "Eta Track; Entries"; numberOfBins = 1000; lowerBound = -3; upperBound = 3;}
823 
824  name_Histogram += add_name;
825  TH1F* histogram = new TH1F(name_Histogram.Data(),discription_Histogram.Data(),numberOfBins,lowerBound,upperBound);
826  histogram->Sumw2();
827  if(j%2==0) histogram->SetLineColor(6);
828  if(j%2==1) histogram->SetLineColor(4);
829  histogram->SetMarkerStyle(20);
830  histogram->SetMarkerSize(0.6);
831  if(j%2==0) histogram->SetMarkerColor(6);
832  if(j%2==1) histogram->SetMarkerColor(4);
833  TH1F* histogram_Clone = (TH1F*)histogram->Clone();
834  listout->Add(histogram_Clone);
835  fDaughterHistogramArray[i][j][k] = histogram_Clone;
836  }
837 
838  TString numberofparticlesperevent="numberofparticlesperevent";
839  numberofparticlesperevent += add_name;
840  TH1F* hist_numberofparticlesperevent = new TH1F(numberofparticlesperevent.Data(),"Number of particles per event; number of particles in one event; Entries",100,0,100);
841  hist_numberofparticlesperevent->Sumw2();
842  hist_numberofparticlesperevent->SetLineColor(6);
843  hist_numberofparticlesperevent->SetMarkerStyle(20);
844  hist_numberofparticlesperevent->SetMarkerSize(0.6);
845  hist_numberofparticlesperevent->SetMarkerColor(6);
846  TH1F* histogram_numberofparticlesperevent = (TH1F*)hist_numberofparticlesperevent->Clone();
847  listout->Add(histogram_numberofparticlesperevent);
848  fDaughterHistogramArray[i][j][10] = histogram_numberofparticlesperevent;
849  }
850 
851  TH1F * effectOfCuts = new TH1F("effectOfCuts","Removal counter",18,0,18);
852  effectOfCuts->SetStats(kTRUE);
853  effectOfCuts->GetXaxis()->SetTitle("Cut number");
854  effectOfCuts->GetYaxis()->SetTitle("Particles cut");
855  effectOfCuts->GetXaxis()->SetBinLabel(1,"total");
856  effectOfCuts->GetXaxis()->SetBinLabel(2,"1");
857  effectOfCuts->GetXaxis()->SetBinLabel(3,"2");
858  effectOfCuts->GetXaxis()->SetBinLabel(4,"3");
859  effectOfCuts->GetXaxis()->SetBinLabel(5,"4");
860  effectOfCuts->GetXaxis()->SetBinLabel(6,"5");
861  effectOfCuts->GetXaxis()->SetBinLabel(7,"6");
862  effectOfCuts->GetXaxis()->SetBinLabel(8,"7");
863  effectOfCuts->GetXaxis()->SetBinLabel(9,"8");
864  effectOfCuts->GetXaxis()->SetBinLabel(10,"9");
865  effectOfCuts->GetXaxis()->SetBinLabel(11,"10");
866  effectOfCuts->GetXaxis()->SetBinLabel(12,"11");
867  effectOfCuts->GetXaxis()->SetBinLabel(13,"12");
868  effectOfCuts->GetXaxis()->SetBinLabel(14,"13");
869  effectOfCuts->GetXaxis()->SetBinLabel(15,"14");
870  effectOfCuts->GetXaxis()->SetBinLabel(16,"15");
871  effectOfCuts->GetXaxis()->SetBinLabel(17,"16");
872  effectOfCuts->GetXaxis()->SetBinLabel(18,"17");
873  listout->Add(effectOfCuts);
874  fDaughterHistogramArrayExtra[i][0] = effectOfCuts;
875 
876  TH1F * effectOfCutsMC = new TH1F("effectOfCutsMC","Removal counter",18,0,18);
877  effectOfCutsMC->SetStats(kTRUE);
878  effectOfCutsMC->GetXaxis()->SetTitle("Cut number");
879  effectOfCutsMC->GetYaxis()->SetTitle("Particles cut");
880  effectOfCutsMC->GetXaxis()->SetBinLabel(1,"total");
881  effectOfCutsMC->GetXaxis()->SetBinLabel(2,"1");
882  effectOfCutsMC->GetXaxis()->SetBinLabel(3,"2");
883  effectOfCutsMC->GetXaxis()->SetBinLabel(4,"3");
884  effectOfCutsMC->GetXaxis()->SetBinLabel(5,"4");
885  effectOfCutsMC->GetXaxis()->SetBinLabel(6,"5");
886  effectOfCutsMC->GetXaxis()->SetBinLabel(7,"6");
887  effectOfCutsMC->GetXaxis()->SetBinLabel(8,"7");
888  effectOfCutsMC->GetXaxis()->SetBinLabel(9,"8");
889  effectOfCutsMC->GetXaxis()->SetBinLabel(10,"9");
890  effectOfCutsMC->GetXaxis()->SetBinLabel(11,"10");
891  effectOfCutsMC->GetXaxis()->SetBinLabel(12,"11");
892  effectOfCutsMC->GetXaxis()->SetBinLabel(13,"12");
893  effectOfCutsMC->GetXaxis()->SetBinLabel(14,"13");
894  effectOfCutsMC->GetXaxis()->SetBinLabel(15,"14");
895  effectOfCutsMC->GetXaxis()->SetBinLabel(16,"15");
896  effectOfCutsMC->GetXaxis()->SetBinLabel(17,"16");
897  effectOfCutsMC->GetXaxis()->SetBinLabel(18,"17");
898  listout->Add(effectOfCutsMC);
899  fDaughterHistogramArrayExtra[i][1] = effectOfCutsMC;
900 
901  TString name_particle_pdg ="particle_pdg";
902  TH1F* hist_particle_pdg = new TH1F(name_particle_pdg.Data(),"Pdg code particle; pdg code; Entries",2000,-0.5,1999.5);
903  hist_particle_pdg->Sumw2();
904  hist_particle_pdg->SetLineColor(6);
905  hist_particle_pdg->SetMarkerStyle(20);
906  hist_particle_pdg->SetMarkerSize(0.6);
907  hist_particle_pdg->SetMarkerColor(6);
908  TH1F* histogram_particle_pdg = (TH1F*)hist_particle_pdg->Clone();
909  listout->Add(histogram_particle_pdg);
910  fDaughterHistogramArrayExtra[i][2] = histogram_particle_pdg;
911 
912  TString name_particle_mother_pdg ="particle_mother_pdg";
913  TH1F* hist_particle_mother_pdg = new TH1F(name_particle_mother_pdg.Data(),"Pdg code particle mother; pdg code; Entries",2000,-0.5,1999.5);
914  hist_particle_mother_pdg->Sumw2();
915  hist_particle_mother_pdg->SetLineColor(6);
916  hist_particle_mother_pdg->SetMarkerStyle(20);
917  hist_particle_mother_pdg->SetMarkerSize(0.6);
918  hist_particle_mother_pdg->SetMarkerColor(6);
919  TH1F* histogram_particle_mother_pdg = (TH1F*)hist_particle_mother_pdg->Clone();
920  listout->Add(histogram_particle_mother_pdg);
921  fDaughterHistogramArrayExtra[i][3] = histogram_particle_mother_pdg;
922 
923  TString name_ptBPlus_vs_ptTrack ="ptBPlus_vs_ptTrackBackground";
924  TH2F* hist_ptBPlus_vs_ptTrack = new TH2F(name_ptBPlus_vs_ptTrack.Data(),"Pt BPlus vs Pt ; p_{T} BPlus [GeV/c]; p_{T} track [GeV/c]",100,0,30,100,0,30);
925  hist_ptBPlus_vs_ptTrack->Sumw2();
926  hist_ptBPlus_vs_ptTrack->SetLineColor(6);
927  hist_ptBPlus_vs_ptTrack->SetMarkerStyle(20);
928  hist_ptBPlus_vs_ptTrack->SetMarkerSize(0.6);
929  hist_ptBPlus_vs_ptTrack->SetMarkerColor(6);
930  TH2F* histogram_ptBPlus_vs_ptTrack = (TH2F*)hist_ptBPlus_vs_ptTrack->Clone();
931  listout->Add(histogram_ptBPlus_vs_ptTrack);
932  fDaughterHistogramArray2D[i][4] = histogram_ptBPlus_vs_ptTrack;
933 
934  TString name_ptBPlus_vs_ptTrackSignal ="ptBPlus_vs_ptTrackSignal";
935  TH2F* hist_ptBPlus_vs_ptTrackSignal = new TH2F(name_ptBPlus_vs_ptTrackSignal.Data(),"Pt BPlus vs Pt ; p_{T} BPlus [GeV/c]; p_{T} track [GeV/c]",100,0,30,100,0,30);
936  hist_ptBPlus_vs_ptTrackSignal->Sumw2();
937  hist_ptBPlus_vs_ptTrackSignal->SetLineColor(4);
938  hist_ptBPlus_vs_ptTrackSignal->SetMarkerStyle(20);
939  hist_ptBPlus_vs_ptTrackSignal->SetMarkerSize(0.6);
940  hist_ptBPlus_vs_ptTrackSignal->SetMarkerColor(6);
941  TH2F* histogram_ptBPlus_vs_ptTrackSignal = (TH2F*)hist_ptBPlus_vs_ptTrackSignal->Clone();
942  listout->Add(histogram_ptBPlus_vs_ptTrackSignal);
943  fDaughterHistogramArray2D[i][5] = histogram_ptBPlus_vs_ptTrackSignal;
944  }
945 
946  //we make the histograms for the reconstructed particles
947  for (Int_t i = 0; i < 3; i++){
948 
949  TString add_name = "";
950  TList * listout;
951  Int_t nHistogramSets = 0;
952  if(i==0) {listout = fOutputD0; nHistogramSets = 6 + 2*fnPtBins;}
953  if(i==1) {listout = fOutputBPlus; nHistogramSets = 6 + 2*fnPtBins;}
954  if(i==2) {listout = fOutputD0_D0Pt; nHistogramSets = 2*fnPtBinsD0forD0ptbin;}
955 
956 
957  for (Int_t j = 0; j < nHistogramSets; j++){
958  if(i<2)
959  {
960  if(j==0) add_name = "";
961  if(j==1) add_name = "Signal";
962  if(j==2) add_name = "Cut";
963  if(j==3) add_name = "SignalCut";
964  if(j==4) add_name = "Result";
965  if(j==5) add_name = "SignalResult";
966  if(j%2==0 && j>5) {add_name = "_ptbin_"; add_name += fPtBinLimits[(j-6)/2]; add_name += "_to_"; add_name += fPtBinLimits[(j-6)/2 + 1];}
967  if(j%2==1 && j>5) {add_name = "Signal_ptbin_"; add_name += fPtBinLimits[(j-7)/2]; add_name += "_to_"; add_name += fPtBinLimits[(j-7)/2 + 1];}
968  }
969  if(i==2)
970  {
971  if(j%2==0) {add_name = "_ptbin_"; add_name += fPtBinLimitsD0forD0ptbin[j/2]; add_name += "_to_"; add_name += fPtBinLimitsD0forD0ptbin[1+j/2];}
972  if(j%2==1) {add_name = "Signal_ptbin_"; add_name += fPtBinLimitsD0forD0ptbin[(j-1)/2]; add_name += "_to_"; add_name += fPtBinLimitsD0forD0ptbin[1+(j-1)/2];}
973  }
974 
975 
976  TString name_Histogram = "";
977  TString discription_Histogram = "";
978  Int_t numberOfBins = 0;
979  Double_t lowerBound = 0.0;
980  Double_t upperBound = 0.0;
981  Int_t numberOfBinsTwo = 0;
982  Double_t lowerBoundTwo = 0.0;
983  Double_t upperBoundTwo = 0.0;
984 
985  for (Int_t k = 0; k < 45; ++k)
986  {
987  if(k==0){name_Histogram = "ptMother"; discription_Histogram = "pt mother; p_{T} [GeV/c]; Entries"; numberOfBins = 300; lowerBound = 0; upperBound = 30;}
988  if(k==1){name_Histogram = "ptFirstDaughter"; discription_Histogram = "pt first daughter; p_{T} [GeV/c]; Entries"; numberOfBins = 300; lowerBound = 0; upperBound = 30;}
989  if(k==2){name_Histogram = "ptSecondDaughter"; discription_Histogram = "pt second daughter; p_{T} [GeV/c]; Entries"; numberOfBins = 300; lowerBound = 0; upperBound = 30;}
990  if(k==3){name_Histogram = "etaMother"; discription_Histogram = "eta mother; #eta; Entries"; numberOfBins = 100; lowerBound = -2; upperBound = 2;}
991  if(k==4){name_Histogram = "phiMother"; discription_Histogram = "phi mother; #phi; Entries"; numberOfBins = 25; lowerBound = 0; upperBound = 2*TMath::Pi();}
992  if(k==5){name_Histogram = "d0Mother"; discription_Histogram = "d0 mother; [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 1.0;}
993  if(k==6){name_Histogram = "d0FirstDaughter"; discription_Histogram = "d0 first daughter; [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 1;}
994 
995  if(k==7){name_Histogram = "d0SecondDaughter"; discription_Histogram = "d0 second daughter; [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 1;}
996 
997  if(k==8){name_Histogram = "pointingAngleMother"; discription_Histogram = "pointing angle; [Cos(#theta)]; Entries"; numberOfBins = 100; lowerBound = -1; upperBound = 1;}
998  if(k==9){name_Histogram = "impactProduct"; discription_Histogram = "impact product; [cm^{2}]; Entries"; numberOfBins = 500; lowerBound = -0.01; upperBound = 0.01;}
999  if(k==10){name_Histogram = "impactProductXY"; discription_Histogram = "impact product XY; [cm^{2}]; Entries"; numberOfBins = 200; lowerBound = 0; upperBound = 0.5;}
1000  if(k==11){name_Histogram = "invariantMassMother"; discription_Histogram = "mass mother candidate; m [GeV/c^{2}]; Entries"; numberOfBins = 10000; lowerBound = 0; upperBound = 10;}
1001  if(k==12){name_Histogram = "deltaMassMother"; discription_Histogram = "mass mother candidate; m [GeV/c^{2}]; Entries"; numberOfBins = 10000; lowerBound = 0; upperBound = 10;}
1002  if(k==13){name_Histogram = "dcaMother"; discription_Histogram = "dca mother; distance [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 0.25;}
1003  if(k==14){name_Histogram = "vertexDistance"; discription_Histogram = "vertex distance between mother and primary vertex; distance [cm]; Entries"; numberOfBins = 100; lowerBound = 0; upperBound = 1;}
1004  if(k==15){name_Histogram = "normDecayLength"; discription_Histogram = "Normalized decay length w.r.t primary vertex; [cm]; Entries"; numberOfBins = 100; lowerBound = 0; upperBound = 50;}
1005  if(k==16){name_Histogram = "pseudoProperDecayTime"; discription_Histogram = "Pseudo Proper Decay Time w.r.t primary vertex; [a.u.]; Entries"; numberOfBins = 1000; lowerBound = -10; upperBound = 10;}
1006  if(k==17){name_Histogram = "DecayTime"; discription_Histogram = "Decay Time w.r.t primary vertex; [a.u.]; Entries"; numberOfBins = 100; lowerBound = 0; upperBound = 0.00000001;}
1007  if(k==18){name_Histogram = "normDecayTime"; discription_Histogram = "Normalized Decay Time w.r.t primary vertex; [a.u.]; Entries"; numberOfBins = 100; lowerBound = 0; upperBound = 0.0000001;}
1008  if(k==19){name_Histogram = "angleMotherFirstDaughter"; discription_Histogram = "flight angle mother and first daughter; [Cos(#phi)]; Entries"; numberOfBins = 100; lowerBound = -1; upperBound = 1;}
1009  if(k==20){name_Histogram = "angleMotherSecondDaughter"; discription_Histogram = "flight angle mother and second daughter; [Cos(#phi)]; Entries"; numberOfBins = 100; lowerBound = 0.5; upperBound = 1;}
1010  if(k==21){name_Histogram = "angleBetweenBothDaughters"; discription_Histogram = "angle between both daughters; [Cos(#phi)]; Entries"; numberOfBins = 100; lowerBound = -1; upperBound = 1;}
1011  if(k==22){name_Histogram = "cosThetaStar"; discription_Histogram = "cosThetaStar; [Cos(#theta*)]; Entries"; numberOfBins = 200; lowerBound = -2; upperBound = 2;}
1012  if(k==23){name_Histogram = "vertexX"; discription_Histogram = "Vertex position; [cm]; Entries"; numberOfBins = 500; lowerBound = -5; upperBound = 5;}
1013  if(k==24){name_Histogram = "vertexY"; discription_Histogram = "Vertex position; [cm]; Entries"; numberOfBins = 500; lowerBound = -5; upperBound = 5;}
1014  if(k==25){name_Histogram = "vertexZ"; discription_Histogram = "Vertex position; [cm]; Entries"; numberOfBins = 500; lowerBound = -20; upperBound = 20;}
1015 
1016 
1017  if(k==26){if(i==0){name_Histogram = "pointingAngleToBPlus"; discription_Histogram = "Pointing angle w.r.t. BPlus decay vertex; [Cos(#theta)]; Entries"; numberOfBins = 200; lowerBound = -1; upperBound = 1;}
1018  else continue;}
1019  if(k==27){if(i==0){name_Histogram = "d0MotherToBPlus"; discription_Histogram = "d0 Mother w.r.t. BPlus decay vertex; [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 1;}
1020  else continue;}
1021  if(k==28){if(i==0){name_Histogram = "d0FirstDaughterToBPlus"; discription_Histogram = "d0 first daughter w.r.t. BPlus decay vertex; [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 1;}
1022  else continue;}
1023  if(k==29){if(i==0){name_Histogram = "d0SecondDaughterToBPlus"; discription_Histogram = "d0 second daughter w.r.t. BPlus decay vertex; [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 1;}
1024  else continue;}
1025  if(k==30){if(i==0){name_Histogram = "impactProductToBPlus"; discription_Histogram = "impact product w.r.t. BPlus decay vertex; [cm]; Entries"; numberOfBins = 400; lowerBound = -0.02; upperBound = 0.02;}
1026  else continue;}
1027  if(k==31){if(i==0){name_Histogram = "impactProductXYToBPlus"; discription_Histogram = "impact product XY w.r.t. BPlus decay vertex; [cm^{2}]; Entries"; numberOfBins = 100; lowerBound = 0; upperBound = 0.5;}
1028  else continue;}
1029  if(k==32){if(i==0){name_Histogram = "normDecayLengthToBPlus"; discription_Histogram = "Normalized decay length w.r.t. BPlus decay vertex; [cm]; Entries"; numberOfBins = 100; lowerBound = 0; upperBound = 50;}
1030  else continue;}
1031  if(k==33){if(i==0){name_Histogram = "pseudoProperDecayTimeToBPlus"; discription_Histogram = "Pseudo Proper Decay Time w.r.t BPlus vertex; [a.u.]; Entries"; numberOfBins = 1000; lowerBound = -1; upperBound = 1;}
1032  else continue;}
1033  if(k==34){if(i==0){name_Histogram = "DecayTimeToBPlus"; discription_Histogram = "Decay Time w.r.t BPlus vertex; [a.u.]; Entries"; numberOfBins = 100; lowerBound = 0; upperBound = 0.00000001;}
1034  else continue;}
1035  if(k==35){if(i==0){name_Histogram = "normDecayTimeToBPlus"; discription_Histogram = "Normalized Decay Time w.r.t BPlus vertex; [a.u.]; Entries"; numberOfBins = 100; lowerBound = 0; upperBound = 0.00000001;}
1036  else continue;}
1037 
1038  if(k==36){name_Histogram = "topomaticFirstDaughter"; discription_Histogram = "topomatic d0 first daughter; [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 20;}
1039  if(k==37){name_Histogram = "topomaticSecondDaughter"; discription_Histogram = "topomatic d0 second daughter; [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 20;}
1040  if(k==38){name_Histogram = "topomaticMax"; discription_Histogram = "Max topomatic; [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 20;}
1041  if(k==39){name_Histogram = "topomaticMin"; discription_Histogram = "Min topomatic; [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 20;}
1042  if(k==40){name_Histogram = "pointingAngleMotherXY"; discription_Histogram = "pointing angle XY; [Cos(#theta)]; Entries"; numberOfBins = 1000; lowerBound = -1; upperBound = 1;}
1043  if(k==41){name_Histogram = "vertexDistanceXY"; discription_Histogram = "vertex distance between mother and primary vertex XY; distance [cm]; Entries"; numberOfBins = 1000; lowerBound = 0; upperBound = 10;}
1044  if(k==42){name_Histogram = "normDecayLengthXY"; discription_Histogram = "Normalized decay length w.r.t primary vertex XY; [cm]; Entries"; numberOfBins = 100; lowerBound = 0; upperBound = 50;}
1045  if(k==43){name_Histogram = "vertexChi2"; discription_Histogram = "#Chi^{2} Vertex; [#Chi^{2}]; Entries"; numberOfBins = 1000; lowerBound = 0; upperBound = 50;}
1046  if(k==44){name_Histogram = "vertexChi2NDF"; discription_Histogram = "#Chi^{2} per NDF Vertex; [#Chi^{2}/NDF]; Entries"; numberOfBins = 1000; lowerBound = 0; upperBound = 50;}
1047 
1048  name_Histogram += add_name;
1049  TH1F* histogram = new TH1F(name_Histogram.Data(),discription_Histogram.Data(),numberOfBins,lowerBound,upperBound);
1050  histogram->Sumw2();
1051  if(j%2==0) histogram->SetLineColor(6);
1052  if(j%2==1) histogram->SetLineColor(4);
1053  histogram->SetMarkerStyle(20);
1054  histogram->SetMarkerSize(0.6);
1055  if(j%2==0) histogram->SetMarkerColor(6);
1056  if(j%2==1) histogram->SetMarkerColor(4);
1057  TH1F* histogram_Clone = (TH1F*)histogram->Clone();
1058  listout->Add(histogram_Clone);
1059  fMotherHistogramArray[i][j][k] = histogram_Clone;
1060  }
1061 
1062 
1063  name_Histogram = "";
1064  discription_Histogram = "";
1065  numberOfBins = 0;
1066  lowerBound = 0.0;
1067  upperBound = 0.0;
1068  numberOfBinsTwo = 0;
1069  lowerBoundTwo = 0.0;
1070  upperBoundTwo = 0.0;
1071 
1072  //we make the 2D histograms for the reconstructed particles
1073  Int_t nFirst = 0;
1074  Int_t nSecond = 1;
1075  Int_t nVariables = 10;
1076  Int_t nHistograms = nVariables * (nVariables - 1) / 2;
1077 
1078  TList * list2D = new TList();
1079  list2D->SetOwner();
1080  TString name2D = "2D_Histograms";
1081  name2D += add_name;
1082  list2D->SetName(name2D);
1083  listout->Add(list2D);
1084 
1085  for (Int_t k = 0; k < nHistograms; ++k)
1086  {
1087  numberOfBins = 50; numberOfBinsTwo = 50;
1088  if(nFirst==0){name_Histogram = "d0FirstDaughter"; discription_Histogram = "d0 first daughter [cm];"; lowerBound = 0; upperBound = 1;}
1089  if(nFirst==1){name_Histogram = "d0SecondDaughter"; discription_Histogram = "d0 second daughter [cm];"; lowerBound = 0; upperBound = 1;}
1090  if(nFirst==2){name_Histogram = "d0Mother"; discription_Histogram = "d0 mother [cm];"; lowerBound = 0; upperBound = 1;}
1091  if(nFirst==3){name_Histogram = "pointingAngleMother"; discription_Histogram = "pointing angle [Cos(#theta)];"; lowerBound = -1; upperBound = 1;}
1092  if(nFirst==4){name_Histogram = "impactProduct"; discription_Histogram = "impact product [cm^{2}];"; lowerBound = -0.01; upperBound = 0.01;}
1093  if(nFirst==5){name_Histogram = "impactProductXY"; discription_Histogram = "impact product XY [cm^{2}];"; lowerBound = 0; upperBound = 0.5;}
1094  if(nFirst==6){name_Histogram = "vertexDistance"; discription_Histogram = "vertex distance between mother and primary vertex [cm];"; lowerBound = 0; upperBound = 1;}
1095  if(nFirst==7){name_Histogram = "normDecayLength"; discription_Histogram = "Normalized decay length w.r.t primary vertex [cm];"; lowerBound = 0; upperBound = 50;}
1096  if(nFirst==8){name_Histogram = "pointingAngleMotherXY"; discription_Histogram = "pointing angle XY [Cos(#theta)];"; lowerBound = -1; upperBound = 1;}
1097  if(nFirst==9){name_Histogram = "vertexDistanceXY"; discription_Histogram = "vertex distance between mother and primary vertex XY [cm];"; lowerBound = 0; upperBound = 1;}
1098  if(nFirst==10){name_Histogram = "normDecayLengthXY"; discription_Histogram = "Normalized decay length w.r.t primary vertex XY [cm];"; lowerBound = 0; upperBound = 50;}
1099 
1100  if(nSecond==0){name_Histogram += "d0FirstDaughter"; discription_Histogram += "d0 first daughter [cm];"; lowerBoundTwo = 0; upperBoundTwo = 1;}
1101  if(nSecond==1){name_Histogram += "d0SecondDaughter"; discription_Histogram += "d0 second daughter [cm];"; lowerBoundTwo = 0; upperBoundTwo = 1;}
1102  if(nSecond==2){name_Histogram += "d0Mother"; discription_Histogram += "d0 mother [cm];"; lowerBoundTwo = 0; upperBoundTwo = 1;}
1103  if(nSecond==3){name_Histogram += "pointingAngleMother"; discription_Histogram += "pointing angle [Cos(#theta)];"; lowerBoundTwo = -1; upperBoundTwo = 1;}
1104  if(nSecond==4){name_Histogram += "impactProduct"; discription_Histogram += "impact product [cm^{2}];"; lowerBoundTwo = -0.01; upperBoundTwo = 0.01;}
1105  if(nSecond==5){name_Histogram += "impactProductXY"; discription_Histogram += "impact product XY [cm^{2}];"; lowerBoundTwo = 0; upperBoundTwo = 0.5;}
1106  if(nSecond==6){name_Histogram += "vertexDistance"; discription_Histogram += "vertex distance between mother and primary vertex [cm];"; lowerBoundTwo = 0; upperBoundTwo = 1;}
1107  if(nSecond==7){name_Histogram += "normDecayLength"; discription_Histogram += "Normalized decay length w.r.t primary vertex [cm];"; lowerBoundTwo = 0; upperBoundTwo = 50;}
1108  if(nSecond==8){name_Histogram += "_pointingAngleMotherXY"; discription_Histogram += "pointing angle XY [Cos(#theta)];"; lowerBoundTwo = -1; upperBoundTwo = 1;}
1109  if(nSecond==9){name_Histogram += "_vertexDistanceXY"; discription_Histogram += "vertex distance between mother and primary vertex XY [cm];"; lowerBoundTwo = 0; upperBoundTwo = 1;}
1110  if(nSecond==10){name_Histogram += "_normDecayLengthXY"; discription_Histogram += "Normalized decay length w.r.t primary vertex XY [cm];"; lowerBoundTwo = 0; upperBoundTwo = 50;}
1111 
1112  name_Histogram += add_name;
1113  TH2F* histogram = new TH2F(name_Histogram.Data(),discription_Histogram.Data(),numberOfBins,lowerBound,upperBound,numberOfBinsTwo,lowerBoundTwo,upperBoundTwo);
1114  histogram->Sumw2();
1115  if(j%2==0) histogram->SetLineColor(6);
1116  if(j%2==1) histogram->SetLineColor(4);
1117  histogram->SetMarkerStyle(20);
1118  histogram->SetMarkerSize(0.6);
1119  histogram->SetMarkerColor(6);
1120  TH2F* histogram_Clone = (TH2F*)histogram->Clone();
1121  list2D->Add(histogram_Clone);
1122  fMotherHistogramArray2D[i][j][k] = histogram_Clone;
1123 
1124  nSecond++;
1125  if(nSecond>nVariables)
1126  {
1127  nFirst++;
1128  nSecond = nFirst + 1;
1129  }
1130  }
1131  }
1132 
1133  TH1F * effectOfCuts = new TH1F("effectOfCuts","Removal counter",100,0,100);
1134  effectOfCuts->SetStats(kTRUE);
1135  effectOfCuts->GetXaxis()->SetTitle("Cut number");
1136  effectOfCuts->GetYaxis()->SetTitle("Particles cut");
1137  effectOfCuts->GetXaxis()->SetBinLabel(1,"total");
1138  for (Int_t i = 1; i < 100; ++i)
1139  {
1140  TString integerText = "";
1141  integerText += i;
1142  effectOfCuts->GetXaxis()->SetBinLabel(i+1,integerText);
1143  }
1144  listout->Add(effectOfCuts);
1145  fMotherHistogramArrayExtra[i][0] = effectOfCuts;
1146 
1147  TH1F * effectOfCutsSignal = new TH1F("effectOfCutsSignal","Removal counter",100,0,100);
1148  effectOfCutsSignal->SetStats(kTRUE);
1149  effectOfCutsSignal->GetXaxis()->SetTitle("Cut number");
1150  effectOfCutsSignal->GetYaxis()->SetTitle("Particles cut");
1151  effectOfCutsSignal->GetXaxis()->SetBinLabel(1,"total");
1152  for (Int_t i = 1; i < 100; ++i)
1153  {
1154  TString integerText = "";
1155  integerText += i;
1156  effectOfCutsSignal->GetXaxis()->SetBinLabel(i+1,integerText);
1157  }
1158  listout->Add(effectOfCutsSignal);
1159  fMotherHistogramArrayExtra[i][1] = effectOfCutsSignal;
1160 
1161  TString name_particle_pdg ="particle_pdg";
1162  TH1F* hist_particle_pdg = new TH1F(name_particle_pdg.Data(),"Pdg code particle; pdg code; Entries",2000,-0.5,1999.5);
1163  hist_particle_pdg->Sumw2();
1164  hist_particle_pdg->SetLineColor(6);
1165  hist_particle_pdg->SetMarkerStyle(20);
1166  hist_particle_pdg->SetMarkerSize(0.6);
1167  hist_particle_pdg->SetMarkerColor(6);
1168  TH1F* histogram_particle_pdg = (TH1F*)hist_particle_pdg->Clone();
1169  listout->Add(histogram_particle_pdg);
1170  fMotherHistogramArrayExtra[i][2] = histogram_particle_pdg;
1171 
1172  TString name_particle_mother_pdg ="particle_mother_pdg";
1173  TH1F* hist_particle_mother_pdg = new TH1F(name_particle_mother_pdg.Data(),"Pdg code particle mother; pdg code; Entries",2000,-0.5,1999.5);
1174  hist_particle_mother_pdg->Sumw2();
1175  hist_particle_mother_pdg->SetLineColor(6);
1176  hist_particle_mother_pdg->SetMarkerStyle(20);
1177  hist_particle_mother_pdg->SetMarkerSize(0.6);
1178  hist_particle_mother_pdg->SetMarkerColor(6);
1179  TH1F* histogram_particle_mother_pdg = (TH1F*)hist_particle_mother_pdg->Clone();
1180  listout->Add(histogram_particle_mother_pdg);
1181  fMotherHistogramArrayExtra[i][3] = histogram_particle_mother_pdg;
1182 
1183  TString name_distance_vertex_from_real ="distance_vertex_from_real";
1184  TH1F* hist_distance_vertex_from_real = new TH1F(name_distance_vertex_from_real.Data(),"Distance reconstructed vertex from real vertex; distance [cm]; Entries",500,0,0.5);
1185  hist_distance_vertex_from_real->Sumw2();
1186  hist_distance_vertex_from_real->SetLineColor(6);
1187  hist_distance_vertex_from_real->SetMarkerStyle(20);
1188  hist_distance_vertex_from_real->SetMarkerSize(0.6);
1189  hist_distance_vertex_from_real->SetMarkerColor(6);
1190  TH1F* histogram_distance_vertex_from_real = (TH1F*)hist_distance_vertex_from_real->Clone();
1191  listout->Add(histogram_distance_vertex_from_real);
1192  fMotherHistogramArrayExtra[i][4] = histogram_distance_vertex_from_real;
1193 
1194  TString name_distance_vertex_from_real_new ="distance_vertex_from_real_new";
1195  TH1F* hist_distance_vertex_from_real_new = new TH1F(name_distance_vertex_from_real_new.Data(),"Distance reconstructed vertex from real vertex; distance [cm]; Entries",500,0,0.5);
1196  hist_distance_vertex_from_real_new->Sumw2();
1197  hist_distance_vertex_from_real_new->SetLineColor(6);
1198  hist_distance_vertex_from_real_new->SetMarkerStyle(20);
1199  hist_distance_vertex_from_real_new->SetMarkerSize(0.6);
1200  hist_distance_vertex_from_real_new->SetMarkerColor(6);
1201  TH1F* histogram_distance_vertex_from_real_new = (TH1F*)hist_distance_vertex_from_real_new->Clone();
1202  listout->Add(histogram_distance_vertex_from_real_new);
1203  fMotherHistogramArrayExtra[i][5] = histogram_distance_vertex_from_real_new;
1204 
1205  TString name_momentum_resolution ="momentum_resolution";
1206  TH1F* hist_momentum_resolution = new TH1F(name_momentum_resolution.Data(),"Momentum resolution; difference between real and reconstructed momentum [GeV/c]; Entries",1000,0,1);
1207  hist_momentum_resolution->Sumw2();
1208  hist_momentum_resolution->SetLineColor(6);
1209  hist_momentum_resolution->SetMarkerStyle(20);
1210  hist_momentum_resolution->SetMarkerSize(0.6);
1211  hist_momentum_resolution->SetMarkerColor(6);
1212  TH1F* histogram_momentum_resolution = (TH1F*)hist_momentum_resolution->Clone();
1213  listout->Add(histogram_momentum_resolution);
1214  fMotherHistogramArrayExtra[i][6] = histogram_momentum_resolution;
1215  }
1216 
1217  //we make the histograms for the same sign method histograms and the pt bins
1218  for (Int_t k = 0; k < fnPtBins+3; ++k){
1219  TString ptBinMother = "";
1220  if(k==0) ptBinMother = "";
1221  if(k==1) ptBinMother = "_ptbin_6_to_inf";
1222  if(k==2) ptBinMother = "_ptbin_3_to_inf";
1223  if(k>2) {ptBinMother += "_ptbin_"; ptBinMother += fPtBinLimits[k-3]; ptBinMother += "_to_"; ptBinMother += fPtBinLimits[k-2];}
1224 
1225  for (Int_t i = 0; i < 8; ++i){
1226  TString signName = "";
1227  if(i==0) signName = "";
1228  if(i==1) signName = "_SameSign";
1229  if(i==2) signName = "_SignSum";
1230  if(i==3) signName = "_HIJING_Background";
1231  if(i==4) signName = "_HIJING_Signal";
1232  if(i==5) signName = "_Background_rotation";
1233  if(i==6) signName = "_HIJING_Background_rotation";
1234  if(i==7) signName = "_correlated511";
1235  TString name_invariantMassMother ="invariantMassBPlus";
1236  name_invariantMassMother += ptBinMother + signName;
1237  TH1F* hist_invariantMassMother = new TH1F(name_invariantMassMother.Data(),"mass mother candidate; m [GeV/c^2]; Entries",2000,0,20);
1238  hist_invariantMassMother->Sumw2();
1239  hist_invariantMassMother->SetLineColor(6);
1240  hist_invariantMassMother->SetMarkerStyle(20);
1241  hist_invariantMassMother->SetMarkerSize(0.6);
1242  hist_invariantMassMother->SetMarkerColor(6);
1243  TH1F* histogram_invariantMassMother = (TH1F*)hist_invariantMassMother->Clone();
1244  fOutputBPlusMC->Add(histogram_invariantMassMother);
1245 
1246  TString name_deltainvariantMassMother ="deltainvariantMassBPlus";
1247  name_deltainvariantMassMother += ptBinMother + signName;
1248  TH1F* hist_deltainvariantMassMother = new TH1F(name_deltainvariantMassMother.Data(),"delta mass mother candidate; m [GeV/c^2]; Entries",2000,0,20);
1249  hist_deltainvariantMassMother->Sumw2();
1250  hist_deltainvariantMassMother->SetLineColor(6);
1251  hist_deltainvariantMassMother->SetMarkerStyle(20);
1252  hist_deltainvariantMassMother->SetMarkerSize(0.6);
1253  hist_deltainvariantMassMother->SetMarkerColor(6);
1254  TH1F* histogram_deltainvariantMassMother = (TH1F*)hist_deltainvariantMassMother->Clone();
1255  fOutputBPlusMC->Add(histogram_deltainvariantMassMother);
1256  }
1257  }
1258 
1259 
1260  TString name_cutEffectBackground ="cutEffectBackground";
1261  TH2I* hist_cutEffectBackground = new TH2I(name_cutEffectBackground.Data(),"Effect of Cuts on background; cut number; cut number",99,0,99,99,0,99);
1262  for (int i = 0; i < 99; ++i)
1263  {
1264  TString integerText = "";
1265  integerText += i;
1266  hist_cutEffectBackground->GetXaxis()->SetBinLabel(i+1,integerText);
1267  hist_cutEffectBackground->GetYaxis()->SetBinLabel(i+1,integerText);
1268  }
1269  TH2I* histogram_cutEffectBackground = (TH2I*)hist_cutEffectBackground->Clone();
1270  fOutputBPlusMC->Add(histogram_cutEffectBackground);
1271 
1272  TString name_cutEffectSignal ="cutEffectSignal";
1273  TH2I* hist_cutEffectSignal = new TH2I(name_cutEffectSignal.Data(),"Effect of Cuts on Signal; cut number; cut number",99,0,99,99,0,99);
1274  for (Int_t i = 0; i < 99; ++i)
1275  {
1276  TString integerText = "";
1277  integerText += i;
1278  hist_cutEffectSignal->GetXaxis()->SetBinLabel(i+1,integerText);
1279  hist_cutEffectSignal->GetYaxis()->SetBinLabel(i+1,integerText);
1280  }
1281  TH2I* histogram_cutEffectSignal = (TH2I*)hist_cutEffectSignal->Clone();
1282  fOutputBPlusMC->Add(histogram_cutEffectSignal);
1283 
1284  TString name_cutEffectUniqueBackground ="cutEffectUniqueBackground";
1285  TH1I* hist_cutEffectUniqueBackground = new TH1I(name_cutEffectUniqueBackground.Data(),"Effect of Cuts on Signal; cut number; cut number",99,0,99);
1286  for (Int_t i = 0; i < 99; ++i)
1287  {
1288  TString integerText = "";
1289  integerText += i;
1290  hist_cutEffectUniqueBackground->GetXaxis()->SetBinLabel(i+1,integerText);
1291  }
1292  TH1I* histogram_cutEffectUniqueBackground = (TH1I*)hist_cutEffectUniqueBackground->Clone();
1293  fOutputBPlusMC->Add(histogram_cutEffectUniqueBackground);
1294 
1295  TString name_cutEffectUniqueSignal ="cutEffectUniqueSignal";
1296  TH1I* hist_cutEffectUniqueSignal = new TH1I(name_cutEffectUniqueSignal.Data(),"Effect of Cuts on Signal; cut number; cut number",99,0,99);
1297  for (Int_t i = 0; i < 99; ++i)
1298  {
1299  TString integerText = "";
1300  integerText += i;
1301  hist_cutEffectUniqueSignal->GetXaxis()->SetBinLabel(i+1,integerText);
1302  }
1303  TH1I* histogram_cutEffectUniqueSignal = (TH1I*)hist_cutEffectUniqueSignal->Clone();
1304  fOutputBPlusMC->Add(histogram_cutEffectUniqueSignal);
1305 
1306  TString name_totalITSBackground ="totalITSBackground";
1307  TH1F* hist_totalITSBackground = new TH1F(name_totalITSBackground.Data(),"Total nr. of ITS hits for the daughters; number [#]; Entries",30,0,30);
1308  hist_totalITSBackground->Sumw2();
1309  hist_totalITSBackground->SetLineColor(6);
1310  hist_totalITSBackground->SetMarkerStyle(20);
1311  hist_totalITSBackground->SetMarkerSize(0.6);
1312  hist_totalITSBackground->SetMarkerColor(6);
1313  TH1F* histogram_totalITSBackground = (TH1F*)hist_totalITSBackground->Clone();
1314  fOutputBPlusMC->Add(histogram_totalITSBackground);
1315 
1316  TString name_totalITSSignal ="totalITSSignal";
1317  TH1F* hist_totalITSSignal = new TH1F(name_totalITSSignal.Data(),"Total nr. of ITS hits for the daughters; number [#]; Entries",30,0,30);
1318  hist_totalITSSignal->Sumw2();
1319  hist_totalITSSignal->SetLineColor(6);
1320  hist_totalITSSignal->SetMarkerStyle(20);
1321  hist_totalITSSignal->SetMarkerSize(0.6);
1322  hist_totalITSSignal->SetMarkerColor(6);
1323  TH1F* histogram_totalITSSignal = (TH1F*)hist_totalITSSignal->Clone();
1324  fOutputBPlusMC->Add(histogram_totalITSSignal);
1325 
1326  TString name_totalTPCBackground ="totalTPCBackground";
1327  TH1F* hist_totalTPCBackground = new TH1F(name_totalTPCBackground.Data(),"Total nr. of TPC hits for the daughters; number [#]; Entries",1000,0,1000);
1328  hist_totalTPCBackground->Sumw2();
1329  hist_totalTPCBackground->SetLineColor(6);
1330  hist_totalTPCBackground->SetMarkerStyle(20);
1331  hist_totalTPCBackground->SetMarkerSize(0.6);
1332  hist_totalTPCBackground->SetMarkerColor(6);
1333  TH1F* histogram_totalTPCBackground = (TH1F*)hist_totalTPCBackground->Clone();
1334  fOutputBPlusMC->Add(histogram_totalTPCBackground);
1335 
1336  TString name_totalTPCSignal ="totalTPCSignal";
1337  TH1F* hist_totalTPCSignal = new TH1F(name_totalTPCSignal.Data(),"Total nr. of TPC hits for the daughters; number [#]; Entries",1000,0,1000);
1338  hist_totalTPCSignal->Sumw2();
1339  hist_totalTPCSignal->SetLineColor(6);
1340  hist_totalTPCSignal->SetMarkerStyle(20);
1341  hist_totalTPCSignal->SetMarkerSize(0.6);
1342  hist_totalTPCSignal->SetMarkerColor(6);
1343  TH1F* histogram_totalTPCSignal = (TH1F*)hist_totalTPCSignal->Clone();
1344  fOutputBPlusMC->Add(histogram_totalTPCSignal);
1345 
1346  TString name_totalSigmaPIDBackground ="totalSigmaPIDBackground";
1347  TH1F* hist_totalSigmaPIDBackground = new TH1F(name_totalSigmaPIDBackground.Data(),"Total sigma of TPC and TOF PID for the daughters; number [#]; Entries",1000,0,100);
1348  hist_totalSigmaPIDBackground->Sumw2();
1349  hist_totalSigmaPIDBackground->SetLineColor(6);
1350  hist_totalSigmaPIDBackground->SetMarkerStyle(20);
1351  hist_totalSigmaPIDBackground->SetMarkerSize(0.6);
1352  hist_totalSigmaPIDBackground->SetMarkerColor(6);
1353  TH1F* histogram_totalSigmaPIDBackground = (TH1F*)hist_totalSigmaPIDBackground->Clone();
1354  fOutputBPlusMC->Add(histogram_totalSigmaPIDBackground);
1355 
1356  TString name_totalSigmaPIDSignal ="totalSigmaPIDSignal";
1357  TH1F* hist_totalSigmaPIDSignal = new TH1F(name_totalSigmaPIDSignal.Data(),"Total sigma of TPC and TOF PID for the daughters; number [#]; Entries",1000,0,100);
1358  hist_totalSigmaPIDSignal->Sumw2();
1359  hist_totalSigmaPIDSignal->SetLineColor(6);
1360  hist_totalSigmaPIDSignal->SetMarkerStyle(20);
1361  hist_totalSigmaPIDSignal->SetMarkerSize(0.6);
1362  hist_totalSigmaPIDSignal->SetMarkerColor(6);
1363  TH1F* histogram_totalSigmaPIDSignal = (TH1F*)hist_totalSigmaPIDSignal->Clone();
1364  fOutputBPlusMC->Add(histogram_totalSigmaPIDSignal);
1365 
1366  TString name_BPlusLabelNumber ="BPlusLabelNumber";
1367  TH1I* hist_BPlusLabelNumber = new TH1I(name_BPlusLabelNumber.Data(),"Number of BPlus mesons with same label; numer with same label; entries",50,0,50);
1368  for (Int_t i = 0; i < 50; ++i)
1369  {
1370  TString integerText = "";
1371  integerText += i;
1372  hist_BPlusLabelNumber->GetXaxis()->SetBinLabel(i+1,integerText);
1373  }
1374  TH1I* histogram_BPlusLabelNumber = (TH1I*)hist_BPlusLabelNumber->Clone();
1375  fOutputBPlusMC->Add(histogram_BPlusLabelNumber);
1376 
1377 
1378  TString name_particle_pdgBPlusPion ="particle_pdgBPlusPion";
1379  TH1F* hist_particle_pdgBPlusPion = new TH1F(name_particle_pdgBPlusPion.Data(),"Pdg code particle; pdg code; Entries",5000,-0.5,4999.5);
1380  hist_particle_pdgBPlusPion->Sumw2();
1381  hist_particle_pdgBPlusPion->SetLineColor(6);
1382  hist_particle_pdgBPlusPion->SetMarkerStyle(20);
1383  hist_particle_pdgBPlusPion->SetMarkerSize(0.6);
1384  hist_particle_pdgBPlusPion->SetMarkerColor(6);
1385  TH1F* histogram_particle_pdgBPlusPion = (TH1F*)hist_particle_pdgBPlusPion->Clone();
1386  fOutputBPlusMC->Add(histogram_particle_pdgBPlusPion);
1387 
1388  TString name_particle_pdgD0First ="particle_pdgD0First";
1389  TH1F* hist_particle_pdgD0First = new TH1F(name_particle_pdgD0First.Data(),"Pdg code particle; pdg code; Entries",5000,-0.5,4999.5);
1390  hist_particle_pdgD0First->Sumw2();
1391  hist_particle_pdgD0First->SetLineColor(6);
1392  hist_particle_pdgD0First->SetMarkerStyle(20);
1393  hist_particle_pdgD0First->SetMarkerSize(0.6);
1394  hist_particle_pdgD0First->SetMarkerColor(6);
1395  TH1F* histogram_particle_pdgD0First = (TH1F*)hist_particle_pdgD0First->Clone();
1396  fOutputBPlusMC->Add(histogram_particle_pdgD0First);
1397 
1398  TString name_particle_pdgD0Second ="particle_pdgD0Second";
1399  TH1F* hist_particle_pdgD0Second = new TH1F(name_particle_pdgD0Second.Data(),"Pdg code particle; pdg code; Entries",5000,-0.5,4999.5);
1400  hist_particle_pdgD0Second->Sumw2();
1401  hist_particle_pdgD0Second->SetLineColor(6);
1402  hist_particle_pdgD0Second->SetMarkerStyle(20);
1403  hist_particle_pdgD0Second->SetMarkerSize(0.6);
1404  hist_particle_pdgD0Second->SetMarkerColor(6);
1405  TH1F* histogram_particle_pdgD0Second = (TH1F*)hist_particle_pdgD0Second->Clone();
1406  fOutputBPlusMC->Add(histogram_particle_pdgD0Second);
1407 
1408  TString name_particle_pdgAll ="particle_pdgAll";
1409  TH1F* hist_particle_pdgAll = new TH1F(name_particle_pdgAll.Data(),"Pdg code particle; pdg code; Entries",5000,-0.5,4999.5);
1410  hist_particle_pdgAll->Sumw2();
1411  hist_particle_pdgAll->SetLineColor(6);
1412  hist_particle_pdgAll->SetMarkerStyle(20);
1413  hist_particle_pdgAll->SetMarkerSize(0.6);
1414  hist_particle_pdgAll->SetMarkerColor(6);
1415  TH1F* histogram_particle_pdgAll = (TH1F*)hist_particle_pdgAll->Clone();
1416  fOutputBPlusMC->Add(histogram_particle_pdgAll);
1417 
1418  TString name_particle_pdgAllSecond ="particle_pdgAllSecond";
1419  TH1F* hist_particle_pdgAllSecond = new TH1F(name_particle_pdgAllSecond.Data(),"Pdg code particle; pdg code; Entries",5000,-0.5,4999.5);
1420  hist_particle_pdgAllSecond->Sumw2();
1421  hist_particle_pdgAllSecond->SetLineColor(6);
1422  hist_particle_pdgAllSecond->SetMarkerStyle(20);
1423  hist_particle_pdgAllSecond->SetMarkerSize(0.6);
1424  hist_particle_pdgAllSecond->SetMarkerColor(6);
1425  TH1F* histogram_particle_pdgAllSecond = (TH1F*)hist_particle_pdgAllSecond->Clone();
1426  fOutputBPlusMC->Add(histogram_particle_pdgAllSecond);
1427 
1428  TString name_particle_pdgAllInvMass ="particle_pdgAllInvMass";
1429  TH2F* hist_particle_pdgAllInvMass = new TH2F(name_particle_pdgAllInvMass.Data(),"Pdg code particle; pdg code; BPlus Inv. Mass",5000,-0.5,4999.5,500,4.0,6.0);
1430  hist_particle_pdgAllInvMass->Sumw2();
1431  hist_particle_pdgAllInvMass->SetLineColor(6);
1432  hist_particle_pdgAllInvMass->SetMarkerStyle(20);
1433  hist_particle_pdgAllInvMass->SetMarkerSize(0.6);
1434  hist_particle_pdgAllInvMass->SetMarkerColor(6);
1435  TH2F* histogram_particle_pdgAllInvMass = (TH2F*)hist_particle_pdgAllInvMass->Clone();
1436  fOutputBPlusMC->Add(histogram_particle_pdgAllInvMass);
1437 
1438  TString name_particle_pdgAllSecondInvMass ="particle_pdgAllSecondInvMass";
1439  TH2F* hist_particle_pdgAllSecondInvMass = new TH2F(name_particle_pdgAllSecondInvMass.Data(),"Pdg code particle; pdg code; BPlus Inv. Mass",5000,-0.5,4999.5,500,4.0,6.0);
1440  hist_particle_pdgAllSecondInvMass->Sumw2();
1441  hist_particle_pdgAllSecondInvMass->SetLineColor(6);
1442  hist_particle_pdgAllSecondInvMass->SetMarkerStyle(20);
1443  hist_particle_pdgAllSecondInvMass->SetMarkerSize(0.6);
1444  hist_particle_pdgAllSecondInvMass->SetMarkerColor(6);
1445  TH2F* histogram_particle_pdgAllSecondInvMass = (TH2F*)hist_particle_pdgAllSecondInvMass->Clone();
1446  fOutputBPlusMC->Add(histogram_particle_pdgAllSecondInvMass);
1447 
1448  TString name_particle_daughterPdgOneStep511 ="particle_daughterPdgOneStep511";
1449  TH2F* hist_particle_daughterPdgOneStep511 = new TH2F(name_particle_daughterPdgOneStep511.Data(),"Pdg daughters; n daughter; Pdg daughter",50,-0.5,49.5,5000,-0.5,4999.5);
1450  hist_particle_daughterPdgOneStep511->Sumw2();
1451  hist_particle_daughterPdgOneStep511->SetLineColor(6);
1452  hist_particle_daughterPdgOneStep511->SetMarkerStyle(20);
1453  hist_particle_daughterPdgOneStep511->SetMarkerSize(0.6);
1454  hist_particle_daughterPdgOneStep511->SetMarkerColor(6);
1455  TH2F* histogram_particle_daughterPdgOneStep511 = (TH2F*)hist_particle_daughterPdgOneStep511->Clone();
1456  fOutputBPlusMC->Add(histogram_particle_daughterPdgOneStep511);
1457 
1458  TString name_particle_daughterPdgOneStep521 ="particle_daughterPdgOneStep521";
1459  TH2F* hist_particle_daughterPdgOneStep521 = new TH2F(name_particle_daughterPdgOneStep521.Data(),"Pdg daughters; n daughter; Pdg daughter",50,-0.5,49.5,5000,-0.5,4999.5);
1460  hist_particle_daughterPdgOneStep521->Sumw2();
1461  hist_particle_daughterPdgOneStep521->SetLineColor(6);
1462  hist_particle_daughterPdgOneStep521->SetMarkerStyle(20);
1463  hist_particle_daughterPdgOneStep521->SetMarkerSize(0.6);
1464  hist_particle_daughterPdgOneStep521->SetMarkerColor(6);
1465  TH2F* histogram_particle_daughterPdgOneStep521 = (TH2F*)hist_particle_daughterPdgOneStep521->Clone();
1466  fOutputBPlusMC->Add(histogram_particle_daughterPdgOneStep521);
1467 
1468  TString name_particle_daughterPdgTwoStep511 ="particle_daughterPdgTwoStep511";
1469  TH2F* hist_particle_daughterPdgTwoStep511 = new TH2F(name_particle_daughterPdgTwoStep511.Data(),"Pdg daughters; n daughter; Pdg daughter",50,-0.5,49.5,5000,-0.5,4999.5);
1470  hist_particle_daughterPdgTwoStep511->Sumw2();
1471  hist_particle_daughterPdgTwoStep511->SetLineColor(6);
1472  hist_particle_daughterPdgTwoStep511->SetMarkerStyle(20);
1473  hist_particle_daughterPdgTwoStep511->SetMarkerSize(0.6);
1474  hist_particle_daughterPdgTwoStep511->SetMarkerColor(6);
1475  TH2F* histogram_particle_daughterPdgTwoStep511 = (TH2F*)hist_particle_daughterPdgTwoStep511->Clone();
1476  fOutputBPlusMC->Add(histogram_particle_daughterPdgTwoStep511);
1477 
1478  TString name_particle_daughterPdgTwoStep521 ="particle_daughterPdgTwoStep521";
1479  TH2F* hist_particle_daughterPdgTwoStep521 = new TH2F(name_particle_daughterPdgTwoStep521.Data(),"Pdg daughters; n daughter; Pdg daughter",50,-0.5,49.5,5000,-0.5,4999.5);
1480  hist_particle_daughterPdgTwoStep521->Sumw2();
1481  hist_particle_daughterPdgTwoStep521->SetLineColor(6);
1482  hist_particle_daughterPdgTwoStep521->SetMarkerStyle(20);
1483  hist_particle_daughterPdgTwoStep521->SetMarkerSize(0.6);
1484  hist_particle_daughterPdgTwoStep521->SetMarkerColor(6);
1485  TH2F* histogram_particle_daughterPdgTwoStep521 = (TH2F*)hist_particle_daughterPdgTwoStep521->Clone();
1486  fOutputBPlusMC->Add(histogram_particle_daughterPdgTwoStep521);
1487 
1488  for (int i = 0; i < 3; ++i)
1489  {
1490  TString name_Histogram;
1491  TString discription_Histogram;
1492  if(i==0) {name_Histogram = "invmassD0PionBPlusPion"; discription_Histogram = "Invariant mass D0 Pion and BPlus Pion; inv. mass [GeV/c^{2}]; Entries";}
1493  if(i==1) {name_Histogram = "invmassD0KaonBPlusPion"; discription_Histogram = "Invariant mass D0 Kaon and BPlus Pion; inv. mass [GeV/c^{2}]; Entries";}
1494  if(i==2) {name_Histogram = "invmassD0PionD0KaonBPlusPion"; discription_Histogram = "Invariant mass D0 Pion, D0 Kaon, and BPlus Pion; inv. mass [GeV/c^{2}]; Entries";}
1495 
1496  for (int j = 0; j < 2; ++j)
1497  {
1498  TString add_name = "";
1499  if(j==1) add_name = "_MC";
1500  name_Histogram += add_name;
1501  TH1F* histogram = new TH1F(name_Histogram.Data(),discription_Histogram.Data(),1000,0,30);
1502  histogram->Sumw2();
1503  if(j%2==0) histogram->SetLineColor(6);
1504  if(j%2==1) histogram->SetLineColor(4);
1505  histogram->SetMarkerStyle(20);
1506  histogram->SetMarkerSize(0.6);
1507  if(j%2==0) histogram->SetMarkerColor(6);
1508  if(j%2==1) histogram->SetMarkerColor(4);
1509  TH1F* histogram_Clone = (TH1F*)histogram->Clone();
1510  fOutputBPlusMC->Add(histogram_Clone);
1511  }
1512  }
1513 
1514 
1515 
1516  return;
1517 }
1518 //-------------------------------------------------------------------------------------
1519 AliAODVertex* AliAnalysisTaskSEBPlustoD0Pi::RecalculateVertex(const AliVVertex *primary,TObjArray *tracks,Double_t bField, Double_t dispersion){
1520  //
1521  // Helper function to recalculate a vertex.
1522  //
1523 
1524  AliESDVertex *vertexESD = 0;
1525  AliAODVertex *vertexAOD = 0;
1526 
1527  AliVertexerTracks vertexer;
1528  vertexer.SetFieldkG(bField);
1529 
1530  vertexer.SetVtxStart((AliESDVertex*)primary); //primary vertex
1531  vertexESD = (AliESDVertex*)vertexer.VertexForSelectedESDTracks(tracks);
1532 
1533  // delete vertexer; vertexer=NULL;
1534 
1535  if(!vertexESD) return vertexAOD;
1536 
1537 
1538  if(vertexESD->GetNContributors()!=tracks->GetEntriesFast())
1539  {
1540  delete vertexESD; vertexESD=nullptr;
1541  return vertexAOD;
1542  }
1543 
1544  // convert to AliAODVertex
1545  Double_t pos[3],cov[6],chi2perNDF;
1546  for(Int_t a=0;a<3;a++)pos[a]=0.;
1547  for(Int_t b=0;b<6;b++)cov[b]=0.;
1548  chi2perNDF=0;
1549 
1550  vertexESD->GetXYZ(pos); // position
1551  vertexESD->GetCovMatrix(cov); //covariance matrix
1552 
1553 
1554  Double_t vertRadius2=pos[0]*pos[0]+pos[1]*pos[1];
1555  if(vertRadius2>8.) //(2.82)^2 radius beam pipe
1556  {
1557  delete vertexESD; vertexESD=nullptr;
1558  return vertexAOD;
1559  }
1560 
1561  chi2perNDF = vertexESD->GetChi2toNDF();
1562  dispersion = vertexESD->GetDispersion();
1563  delete vertexESD; vertexESD=nullptr;
1564  Int_t nprongs = 2; //tracks->GetEntriesFast();
1565  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
1566 
1567  return vertexAOD;
1568 }
1569 //-------------------------------------------------------------------------------------
1570 void AliAnalysisTaskSEBPlustoD0Pi::BPlustoD0PiSignalTracksInMC(TClonesArray * mcTrackArray,AliAODEvent* aodevent,TMatrix * BPlustoD0PiLabelMatrix, TList *listout){
1571 
1572  TMatrix &particleMatrix = *BPlustoD0PiLabelMatrix;
1573  for (Int_t i=0; i<mcTrackArray->GetEntriesFast(); i++){
1574 
1575  Int_t mcLabelPionBPlus = 0;
1576  Int_t mcLabelPionD0 = 0;
1577  Int_t mcLabelKaon = 0;
1578  Int_t mcLabelD0 = 0;
1579  Int_t mcLabelBPlus = 0;
1580 
1581  Double_t ptMC[5] = {0.0};
1582  Double_t yMC[5] = {0.0};
1583  Double_t pseudoYMC[5] = {0.0};
1584 
1585  Bool_t mcPionBPlusPresent = kFALSE;
1586  Bool_t mcPionD0Present = kFALSE;
1587  Bool_t mcKaonPresent = kFALSE;
1588 
1589 
1590  AliAODMCParticle *mcTrackParticle = dynamic_cast< AliAODMCParticle*>(mcTrackArray->At(i));
1591  if(!mcTrackParticle) {std::cout << "no particle" << std::endl; continue;}
1592  Int_t pdgCodeMC=TMath::Abs(mcTrackParticle->GetPdgCode());
1593 
1594  if (pdgCodeMC==521)
1595  { //if the track is a BPlus we look at its daughters
1596 
1597  mcLabelBPlus = i;
1598  Int_t nDaughterBPlus = mcTrackParticle->GetNDaughters();
1599  ptMC[0] = mcTrackParticle->Pt();
1600  yMC[0] = mcTrackParticle->Y();
1601  pseudoYMC[0] = mcTrackParticle->Eta();
1602 
1603  TString fillthis= "BPlus_in_analysis";
1604  ((TH1F*)(listout->FindObject(fillthis)))->Fill(0);
1605 
1606  if(nDaughterBPlus==2)
1607  {
1608  for(Int_t iDaughterBPlus=0; iDaughterBPlus<2; iDaughterBPlus++)
1609  {
1610  AliAODMCParticle* daughterBPlus = (AliAODMCParticle*)mcTrackArray->At(mcTrackParticle->GetDaughter(iDaughterBPlus));
1611  if(!daughterBPlus) break;
1612  Int_t pdgCodeDaughterBPlus=TMath::Abs(daughterBPlus->GetPdgCode());
1613 
1614  if (pdgCodeDaughterBPlus==211) //if the track is a pion we save its monte carlo label
1615  {
1616  mcLabelPionBPlus = mcTrackParticle->GetDaughter(iDaughterBPlus);
1617  mcPionBPlusPresent = kTRUE;
1618  ptMC[1] = daughterBPlus->Pt();
1619  yMC[1] = daughterBPlus->Y();
1620  pseudoYMC[1] = daughterBPlus->Eta();
1621  } else if (pdgCodeDaughterBPlus==421) //if the track is a D0 we look at its daughters
1622  {
1623  mcLabelD0 = mcTrackParticle->GetDaughter(iDaughterBPlus);
1624  Int_t nDaughterD0 = daughterBPlus->GetNDaughters();
1625  ptMC[2] = daughterBPlus->Pt();
1626  yMC[2] = daughterBPlus->Y();
1627  pseudoYMC[2] = daughterBPlus->Eta();
1628 
1629  if(nDaughterD0==2)
1630  {
1631  for(Int_t iDaughterD0=0; iDaughterD0<2; iDaughterD0++)
1632  {
1633  AliAODMCParticle* daughterD0 = (AliAODMCParticle*)mcTrackArray->At(daughterBPlus->GetDaughter(iDaughterD0));
1634  if(!daughterD0) break;
1635  Int_t pdgCodeDaughterD0=TMath::Abs(daughterD0->GetPdgCode());
1636  if (pdgCodeDaughterD0==211) //if the track is a pion we save its monte carlo label
1637  {
1638  mcLabelPionD0 = daughterBPlus->GetDaughter(iDaughterD0);
1639  ptMC[3] = daughterD0->Pt();
1640  yMC[3] = daughterD0->Y();
1641  pseudoYMC[3] = daughterD0->Eta();
1642  mcPionD0Present = kTRUE;
1643  } else if (pdgCodeDaughterD0==321) //if the track is a kaon we save its monte carlo label
1644  {
1645  mcLabelKaon = daughterBPlus->GetLabel();
1646  mcKaonPresent = kTRUE;
1647  ptMC[4] = daughterD0->Pt();
1648  yMC[4] = daughterD0->Y();
1649  pseudoYMC[4] = daughterD0->Eta();
1650  } else break;
1651  }
1652  }
1653  } else break;
1654  }
1655  }
1656  }
1657  if (mcPionBPlusPresent && mcPionD0Present && mcKaonPresent){
1658 
1659  TString fillthis= "BPlus_in_analysis";
1660  ((TH1F*)(listout->FindObject(fillthis)))->Fill(1);
1661 
1662  fillthis= "BPlus_per_bin";
1663  for (Int_t j = 0; j < fnPtBins; ++j)
1664  {
1665  if(fPtBinLimits[j] < ptMC[0] && ptMC[0] < fPtBinLimits[j+1]) {((TH1F*)(listout->FindObject(fillthis)))->Fill(j); break;}
1666  }
1667 
1668 
1669  fillthis= "mc_BPlus_pt";
1670  ((TH1F*)(listout->FindObject(fillthis)))->Fill(ptMC[0]);
1671  fillthis= "mc_BPlus_pion_pt";
1672  ((TH1F*)(listout->FindObject(fillthis)))->Fill(ptMC[1]);
1673  fillthis= "mc_D0_pt";
1674  ((TH1F*)(listout->FindObject(fillthis)))->Fill(ptMC[2]);
1675  fillthis= "mc_D0_pion_pt";
1676  ((TH1F*)(listout->FindObject(fillthis)))->Fill(ptMC[3]);
1677  fillthis= "mc_D0_kaon_pt";
1678  ((TH1F*)(listout->FindObject(fillthis)))->Fill(ptMC[4]);
1679 
1680  fillthis= "mc_BPlus_rapidity_true";
1681  ((TH1F*)(listout->FindObject(fillthis)))->Fill(yMC[0]);
1682  fillthis= "mc_BPlus_pion_rapidity_true";
1683  ((TH1F*)(listout->FindObject(fillthis)))->Fill(yMC[1]);
1684  fillthis= "mc_D0_rapidity_true";
1685  ((TH1F*)(listout->FindObject(fillthis)))->Fill(yMC[2]);
1686  fillthis= "mc_D0_pion_rapidity_true";
1687  ((TH1F*)(listout->FindObject(fillthis)))->Fill(yMC[3]);
1688  fillthis= "mc_D0_kaon_rapidity_true";
1689  ((TH1F*)(listout->FindObject(fillthis)))->Fill(yMC[4]);
1690 
1691  fillthis= "mc_BPlus_pseudorapidity_true";
1692  ((TH1F*)(listout->FindObject(fillthis)))->Fill(pseudoYMC[0]);
1693  fillthis= "mc_BPlus_pion_pseudorapidity_true";
1694  ((TH1F*)(listout->FindObject(fillthis)))->Fill(pseudoYMC[1]);
1695  fillthis= "mc_D0_pseudorapidity_true";
1696  ((TH1F*)(listout->FindObject(fillthis)))->Fill(pseudoYMC[2]);
1697  fillthis= "mc_D0_pion_pseudorapidity_true";
1698  ((TH1F*)(listout->FindObject(fillthis)))->Fill(pseudoYMC[3]);
1699  fillthis= "mc_D0_kaon_pseudorapidity_true";
1700  ((TH1F*)(listout->FindObject(fillthis)))->Fill(pseudoYMC[4]);
1701 
1702  // We check if the tracks are in acceptance
1703  if(ptMC[1] < 0.1 || TMath::Abs(pseudoYMC[1]) > 0.8 ) continue;
1704  if(ptMC[3] < 0.1 || TMath::Abs(pseudoYMC[3]) > 0.8 ) continue;
1705  if(ptMC[4] < 0.1 || TMath::Abs(pseudoYMC[4]) > 0.8 ) continue;
1706 
1707  // We check if the BPlus is in the fiducial region
1708  if(TMath::Abs(yMC[0]) > 0.8) continue;
1709 
1710  Int_t rows = BPlustoD0PiLabelMatrix->GetNrows();
1711 
1712  BPlustoD0PiLabelMatrix->ResizeTo(rows+1,5);
1713  particleMatrix(rows,0) = mcLabelPionBPlus;
1714  particleMatrix(rows,1) = mcLabelPionD0;
1715  particleMatrix(rows,2) = mcLabelKaon;
1716  particleMatrix(rows,3) = mcLabelD0;
1717  particleMatrix(rows,4) = mcLabelBPlus;
1718 
1719  fillthis= "BPlus_in_analysis";
1720  ((TH1F*)(listout->FindObject(fillthis)))->Fill(2);
1721 
1722 
1723  fillthis= "BPlus_per_bin_in_Acc";
1724  for (Int_t j = 0; j < fnPtBins; ++j)
1725  {
1726  if(fPtBinLimits[j] < ptMC[0] && ptMC[0] < fPtBinLimits[j+1]) {((TH1F*)(listout->FindObject(fillthis)))->Fill(j); break;}
1727  }
1728  }
1729  }
1730 
1731 
1732  return;
1733 }
1734 //-------------------------------------------------------------------------------------
1735 Bool_t AliAnalysisTaskSEBPlustoD0Pi::D0FirstDaughterSelection(AliAODTrack* aodTrack, AliAODVertex *primaryVertex, Double_t bz, TClonesArray * mcTrackArray, TMatrix * BPlustoD0PiLabelMatrix, AliAODMCHeader * header){
1736 
1737  // we select the D0 pion and save its information
1738  if(!aodTrack) AliFatal("Not a standard AOD");
1739 
1740  //quick quality cut
1741  if(aodTrack->GetITSNcls() < 1) return kFALSE;
1742  if(aodTrack->GetTPCNcls() < 1) return kFALSE;
1743  if(aodTrack->GetStatus()&AliESDtrack::kITSpureSA) return kFALSE;
1744  if(!(aodTrack->GetStatus()&AliESDtrack::kITSin)) return kFALSE;
1745  if(aodTrack->GetID() < 0) return kFALSE;
1746  Double_t covtest[21];
1747  if(!aodTrack->GetCovarianceXYZPxPyPz(covtest)) return kFALSE;
1748 
1749  Int_t mcLabelParticle = -1;
1750  Int_t pdgParticle = -1;
1751  mcLabelParticle = aodTrack->GetLabel();
1752 
1753  // we fill histograms with information of the track
1754  Double_t pt_track = aodTrack->Pt();
1755  Double_t momentum_track = aodTrack->P();
1756  Int_t numberOfITS = aodTrack->GetITSNcls();
1757  Int_t numberOfTPC = aodTrack->GetTPCNcls();
1758 
1759  AliExternalTrackParam particleTrack;
1760  particleTrack.CopyFromVTrack(aodTrack);
1761  Double_t d0[2],covd0[3];
1762  particleTrack.PropagateToDCA(primaryVertex,bz,100.,d0,covd0);
1763 
1764  //we check if the particle is a signal track, we look at both daughter options therefore the signal will be too high in the D0 daughter signal histograms
1765  Bool_t isDesiredCandidate = kFALSE;
1766  if(fUseMCInfo){
1767  TMatrix &particleMatrix = *BPlustoD0PiLabelMatrix;
1768  for (Int_t k = 0; k < BPlustoD0PiLabelMatrix->GetNrows(); ++k){
1769  if(mcLabelParticle == (Int_t)particleMatrix(k,1) || mcLabelParticle == (Int_t)particleMatrix(k,2)){
1770  isDesiredCandidate = kTRUE;
1771  break;
1772  }
1773  }
1774  }
1775 
1776  if(fUseMCInfo){
1777  if(IsTrackInjected(aodTrack,header,mcTrackArray) && !isDesiredCandidate && fQuickSignalAnalysis == 2) return kFALSE;
1778  }
1779 
1780  Int_t daughterType = 0;
1781 
1782 
1783  Int_t histType = 0;
1784  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
1785  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
1786  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
1787  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
1788 
1789  for (Int_t j = 0; j < 10; ++j)
1790  {
1791  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
1792 
1793  }
1794  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
1795 
1796  if(isDesiredCandidate)
1797  {
1798  histType = 1;
1799  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
1800  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
1801  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
1802  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
1803 
1804  for (Int_t j = 0; j < 10; ++j)
1805  {
1806  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
1807 
1808  }
1809  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
1810  }
1811 
1812  //we apply a number of cuts on the particle
1813  Bool_t bCut = kFALSE;
1814 
1815  //We do not apply a PID cut at this stage since we don't know if we are dealing with a kaon or a pion
1816 
1817  if(aodTrack->GetITSNcls() < fCuts->GetMinITSNclsD0FirstDaughter()){
1818  if(isDesiredCandidate) {
1819  ((TH1F*)fDaughterHistogramArrayExtra[0][1])->Fill(3);
1820  } else ((TH1F*)fDaughterHistogramArrayExtra[0][0])->Fill(3);
1821  bCut = kTRUE;
1822  }
1823 
1824  if(aodTrack->GetTPCNcls() < fCuts->GetMinTPCNclsD0FirstDaughter()){
1825  if(isDesiredCandidate) {
1826  ((TH1F*)fDaughterHistogramArrayExtra[0][1])->Fill(4);
1827  } else ((TH1F*)fDaughterHistogramArrayExtra[0][0])->Fill(4);
1828  bCut = kTRUE;
1829  }
1830 
1831  if(fCuts->UseITSRefitD0FirstDaughter()==kTRUE){
1832  if(!(aodTrack->GetStatus()&AliESDtrack::kITSrefit)) {
1833  if(isDesiredCandidate) {
1834  ((TH1F*)fDaughterHistogramArrayExtra[0][1])->Fill(5);
1835  } else ((TH1F*)fDaughterHistogramArrayExtra[0][0])->Fill(5);
1836  bCut = kTRUE;
1837  }
1838  }
1839 
1840  if(fCuts->UseTPCRefitD0FirstDaughter()==kTRUE){
1841  if((!(aodTrack->GetStatus()&AliESDtrack::kTPCrefit))) {
1842  if(isDesiredCandidate) {
1843  ((TH1F*)fDaughterHistogramArrayExtra[0][1])->Fill(6);
1844  } else ((TH1F*)fDaughterHistogramArrayExtra[0][0])->Fill(6);
1845  bCut = kTRUE;
1846  }
1847  }
1848 
1849  if(fCuts->UseFilterBitD0FirstDaughter()==kTRUE){
1850  if(!(aodTrack->TestFilterMask(BIT(fCuts->GetFilterBitD0FirstDaughter())))) {
1851  if(isDesiredCandidate) {
1852  ((TH1F*)fDaughterHistogramArrayExtra[0][1])->Fill(7);
1853  } else ((TH1F*)fDaughterHistogramArrayExtra[0][0])->Fill(7);
1854  bCut = kTRUE;
1855  }
1856  }
1857 
1858  if(aodTrack->Pt() < fCuts->GetMinPtD0FirstDaughter()){
1859  if(isDesiredCandidate) {
1860  ((TH1F*)fDaughterHistogramArrayExtra[0][1])->Fill(8);
1861  } else ((TH1F*)fDaughterHistogramArrayExtra[0][0])->Fill(8);
1862  bCut = kTRUE;
1863  }
1864 
1865  if(TMath::Abs(d0[0]) < fCuts->GetMind0D0FirstDaughter()){
1866  if(isDesiredCandidate) {
1867  ((TH1F*)fDaughterHistogramArrayExtra[0][1])->Fill(12);
1868  } else ((TH1F*)fDaughterHistogramArrayExtra[0][0])->Fill(12);
1869  bCut = kTRUE;
1870  }
1871 
1872  if(TMath::Abs(aodTrack->Eta()) > fCuts->GetMaxAbsEtaD0FirstDaughter()){
1873  if(isDesiredCandidate) {
1874  ((TH1F*)fDaughterHistogramArrayExtra[0][1])->Fill(9);
1875  } else ((TH1F*)fDaughterHistogramArrayExtra[0][0])->Fill(9);
1876  bCut = kTRUE;
1877  }
1878 
1879  Bool_t bHardSelectionArrayITS[7] = {kFALSE};
1880  fCuts->GetHardSelectionArrayITSD0FirstDaughter(bHardSelectionArrayITS);
1881  Bool_t bSoftSelectionArrayITS[7] = {kFALSE};
1882  fCuts->GetSoftSelectionArrayITSD0FirstDaughter(bSoftSelectionArrayITS);
1883 
1884  Bool_t bHardITSPass = kTRUE;
1885  for (Int_t j = 0; j < 7; ++j)
1886  {
1887  if(bHardSelectionArrayITS[j])
1888  {
1889  if(!aodTrack->HasPointOnITSLayer(j)) bHardITSPass = kFALSE;
1890  }
1891  }
1892 
1893  Int_t nCounterSoftSelection = 0;
1894  Bool_t bSoftITSPass = kTRUE;
1895  for (Int_t j = 0; j < 7; ++j)
1896  {
1897  if(bSoftSelectionArrayITS[j])
1898  {
1899  if(aodTrack->HasPointOnITSLayer(j)) nCounterSoftSelection++;
1900  }
1901  }
1902  if(nCounterSoftSelection < fCuts->GetNSoftITSCutD0FirstDaughter()) bSoftITSPass = kFALSE;
1903 
1904  if(!bHardITSPass){
1905  if(isDesiredCandidate) {
1906  ((TH1F*)fDaughterHistogramArrayExtra[0][1])->Fill(10);
1907  } else ((TH1F*)fDaughterHistogramArrayExtra[0][0])->Fill(10);
1908  bCut = kTRUE;
1909  }
1910 
1911  if(!bSoftITSPass){
1912  if(isDesiredCandidate) {
1913  ((TH1F*)fDaughterHistogramArrayExtra[0][1])->Fill(11);
1914  } else ((TH1F*)fDaughterHistogramArrayExtra[0][0])->Fill(11);
1915  bCut = kTRUE;
1916  }
1917 
1918  if(!isDesiredCandidate && fQuickSignalAnalysis == 1) bCut = kTRUE;
1919 
1920  if(bCut) {
1921  if(isDesiredCandidate) {
1922  ((TH1F*)fDaughterHistogramArrayExtra[0][1])->Fill(0);
1923  } else ((TH1F*)fDaughterHistogramArrayExtra[0][0])->Fill(0);
1924  return kFALSE;
1925  }
1926 
1927  //we fill histograms with track information of the tracks that pass the cuts
1928  histType = 2;
1929  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
1930  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
1931  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
1932  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
1933 
1934  for (Int_t j = 0; j < 10; ++j)
1935  {
1936  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
1937 
1938  }
1939  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
1940 
1941  if(isDesiredCandidate)
1942  {
1943  histType = 3;
1944  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
1945  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
1946  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
1947  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
1948 
1949  for (Int_t j = 0; j < 10; ++j)
1950  {
1951  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
1952 
1953  }
1954  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
1955  }
1956 
1957  return kTRUE;
1958 }
1959 //-------------------------------------------------------------------------------------
1960 Bool_t AliAnalysisTaskSEBPlustoD0Pi::D0SecondDaughterSelection(AliAODTrack* aodTrack, AliAODVertex *primaryVertex, Double_t bz, TClonesArray * mcTrackArray, TMatrix * BPlustoD0PiLabelMatrix, AliAODMCHeader * header){
1961 
1962  // we select the D0 pion and save its information
1963  if(!aodTrack) AliFatal("Not a standard AOD");
1964 
1965  //quick quality cut
1966  if(aodTrack->GetITSNcls() < 1) return kFALSE;
1967  if(aodTrack->GetTPCNcls() < 1) return kFALSE;
1968  if(aodTrack->GetStatus()&AliESDtrack::kITSpureSA) return kFALSE;
1969  if(!(aodTrack->GetStatus()&AliESDtrack::kITSin)) return kFALSE;
1970  if(aodTrack->GetID() < 0) return kFALSE;
1971  Double_t covtest[21];
1972  if(!aodTrack->GetCovarianceXYZPxPyPz(covtest)) return kFALSE;
1973 
1974  Int_t mcLabelParticle = -1;
1975  Int_t pdgParticle = -1;
1976  mcLabelParticle = aodTrack->GetLabel();
1977 
1978  // we fill histograms with information of the track
1979  Double_t pt_track = aodTrack->Pt();
1980  Double_t momentum_track = aodTrack->P();
1981  Int_t numberOfITS = aodTrack->GetITSNcls();
1982  Int_t numberOfTPC = aodTrack->GetTPCNcls();
1983 
1984  AliExternalTrackParam particleTrack;
1985  particleTrack.CopyFromVTrack(aodTrack);
1986  Double_t d0[2],covd0[3];
1987  particleTrack.PropagateToDCA(primaryVertex,bz,100.,d0,covd0);
1988 
1989  //we check if the particle is a signal track, we look at both daughter options therefore the signal will be too high in the D0 daughter signal histograms
1990  Bool_t isDesiredCandidate = kFALSE;
1991  if(fUseMCInfo){
1992  TMatrix &particleMatrix = *BPlustoD0PiLabelMatrix;
1993  for (Int_t k = 0; k < BPlustoD0PiLabelMatrix->GetNrows(); ++k){
1994  if(mcLabelParticle == (Int_t)particleMatrix(k,1) || mcLabelParticle == (Int_t)particleMatrix(k,2)){
1995  isDesiredCandidate = kTRUE;
1996  break;
1997  }
1998  }
1999  }
2000 
2001  if(fUseMCInfo){
2002  if(IsTrackInjected(aodTrack,header,mcTrackArray) && !isDesiredCandidate && fQuickSignalAnalysis == 2) return kFALSE;
2003  }
2004 
2005  Int_t daughterType = 1;
2006 
2007 
2008  Int_t histType = 0;
2009  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
2010  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
2011  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
2012  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
2013 
2014  for (Int_t j = 0; j < 10; ++j)
2015  {
2016  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
2017 
2018  }
2019  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
2020 
2021  if(isDesiredCandidate)
2022  {
2023  histType = 1;
2024  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
2025  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
2026  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
2027  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
2028 
2029  for (Int_t j = 0; j < 10; ++j)
2030  {
2031  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
2032 
2033  }
2034  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
2035  }
2036 
2037  //we apply a number of cuts on the particle
2038  Bool_t bCut = kFALSE;
2039 
2040  //We do not apply a PID cut at this stage since we don't know if we are dealing with a kaon or a pion
2041 
2042  if(aodTrack->GetITSNcls() < fCuts->GetMinITSNclsD0SecondDaughter()){
2043  if(isDesiredCandidate) {
2044  ((TH1F*)fDaughterHistogramArrayExtra[1][1])->Fill(3);
2045  } else ((TH1F*)fDaughterHistogramArrayExtra[1][0])->Fill(3);
2046  bCut = kTRUE;
2047  }
2048 
2049  if(aodTrack->GetTPCNcls() < fCuts->GetMinTPCNclsD0SecondDaughter()){
2050  if(isDesiredCandidate) {
2051  ((TH1F*)fDaughterHistogramArrayExtra[1][1])->Fill(4);
2052  } else ((TH1F*)fDaughterHistogramArrayExtra[1][0])->Fill(4);
2053  bCut = kTRUE;
2054  }
2055 
2056  if(fCuts->UseITSRefitD0SecondDaughter()==kTRUE){
2057  if(!(aodTrack->GetStatus()&AliESDtrack::kITSrefit)) {
2058  if(isDesiredCandidate) {
2059  ((TH1F*)fDaughterHistogramArrayExtra[1][1])->Fill(5);
2060  } else ((TH1F*)fDaughterHistogramArrayExtra[1][0])->Fill(5);
2061  bCut = kTRUE;
2062  }
2063  }
2064 
2065  if(fCuts->UseTPCRefitD0SecondDaughter()==kTRUE){
2066  if((!(aodTrack->GetStatus()&AliESDtrack::kTPCrefit))) {
2067  if(isDesiredCandidate) {
2068  ((TH1F*)fDaughterHistogramArrayExtra[1][1])->Fill(6);
2069  } else ((TH1F*)fDaughterHistogramArrayExtra[1][0])->Fill(6);
2070  bCut = kTRUE;
2071  }
2072  }
2073 
2074  if(fCuts->UseFilterBitD0SecondDaughter()==kTRUE){
2075  if(!(aodTrack->TestFilterMask(BIT(fCuts->GetFilterBitD0SecondDaughter())))) {
2076  if(isDesiredCandidate) {
2077  ((TH1F*)fDaughterHistogramArrayExtra[1][1])->Fill(7);
2078  } else ((TH1F*)fDaughterHistogramArrayExtra[1][0])->Fill(7);
2079  bCut = kTRUE;
2080  }
2081  }
2082 
2083  if(aodTrack->Pt() < fCuts->GetMinPtD0SecondDaughter()){
2084  if(isDesiredCandidate) {
2085  ((TH1F*)fDaughterHistogramArrayExtra[1][1])->Fill(8);
2086  } else ((TH1F*)fDaughterHistogramArrayExtra[1][0])->Fill(8);
2087  bCut = kTRUE;
2088  }
2089 
2090  if(TMath::Abs(d0[0]) < fCuts->GetMind0D0SecondDaughter()){
2091  if(isDesiredCandidate) {
2092  ((TH1F*)fDaughterHistogramArrayExtra[1][1])->Fill(12);
2093  } else ((TH1F*)fDaughterHistogramArrayExtra[1][0])->Fill(12);
2094  bCut = kTRUE;
2095  }
2096 
2097  if(TMath::Abs(aodTrack->Eta()) > fCuts->GetMaxAbsEtaD0SecondDaughter()){
2098  if(isDesiredCandidate) {
2099  ((TH1F*)fDaughterHistogramArrayExtra[1][1])->Fill(9);
2100  } else ((TH1F*)fDaughterHistogramArrayExtra[1][0])->Fill(9);
2101  bCut = kTRUE;
2102  }
2103 
2104  Bool_t bHardSelectionArrayITS[7] = {kFALSE};
2105  fCuts->GetHardSelectionArrayITSD0SecondDaughter(bHardSelectionArrayITS);
2106  Bool_t bSoftSelectionArrayITS[7] = {kFALSE};
2107  fCuts->GetSoftSelectionArrayITSD0SecondDaughter(bSoftSelectionArrayITS);
2108 
2109  Bool_t bHardITSPass = kTRUE;
2110  for (Int_t j = 0; j < 7; ++j)
2111  {
2112  if(bHardSelectionArrayITS[j])
2113  {
2114  if(!aodTrack->HasPointOnITSLayer(j)) bHardITSPass = kFALSE;
2115  }
2116  }
2117 
2118  Int_t nCounterSoftSelection = 0;
2119  Bool_t bSoftITSPass = kTRUE;
2120  for (Int_t j = 0; j < 7; ++j)
2121  {
2122  if(bSoftSelectionArrayITS[j])
2123  {
2124  if(aodTrack->HasPointOnITSLayer(j)) nCounterSoftSelection++;
2125  }
2126  }
2127  if(nCounterSoftSelection < fCuts->GetNSoftITSCutD0SecondDaughter()) bSoftITSPass = kFALSE;
2128 
2129  if(!bHardITSPass){
2130  if(isDesiredCandidate) {
2131  ((TH1F*)fDaughterHistogramArrayExtra[1][1])->Fill(10);
2132  } else ((TH1F*)fDaughterHistogramArrayExtra[1][0])->Fill(10);
2133  bCut = kTRUE;
2134  }
2135 
2136  if(!bSoftITSPass){
2137  if(isDesiredCandidate) {
2138  ((TH1F*)fDaughterHistogramArrayExtra[1][1])->Fill(11);
2139  } else ((TH1F*)fDaughterHistogramArrayExtra[1][0])->Fill(11);
2140  bCut = kTRUE;
2141  }
2142 
2143  if(!isDesiredCandidate && fQuickSignalAnalysis == 1) bCut = kTRUE;
2144 
2145  if(bCut) {
2146  if(isDesiredCandidate) {
2147  ((TH1F*)fDaughterHistogramArrayExtra[1][1])->Fill(0);
2148  } else ((TH1F*)fDaughterHistogramArrayExtra[1][0])->Fill(0);
2149  return kFALSE;
2150  }
2151 
2152  //we fill histograms with track information of the tracks that pass the cuts
2153  histType = 2;
2154  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
2155  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
2156  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
2157  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
2158 
2159  for (Int_t j = 0; j < 10; ++j)
2160  {
2161  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
2162 
2163  }
2164  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
2165 
2166  if(isDesiredCandidate)
2167  {
2168  histType = 3;
2169  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
2170  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
2171  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
2172  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
2173 
2174  for (Int_t j = 0; j < 10; ++j)
2175  {
2176  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
2177 
2178  }
2179  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
2180  }
2181 
2182  return kTRUE;
2183 }
2184 //-------------------------------------------------------------------------------------
2185 void AliAnalysisTaskSEBPlustoD0Pi::BPlusPionSelection(AliAODEvent* aodEvent, AliAODVertex *primaryVertex, Double_t bz, TClonesArray * mcTrackArray, TMatrix * BPlustoD0PiLabelMatrix, AliAODMCHeader * header){
2186 
2187  //we keep track of the number of particles we could use and how many we actually use after cuts
2188  Int_t numberofparticles = 0;
2189  Int_t numberofparticlesused = 0;
2190  Int_t iClonesArray = 0;
2191 
2192  TString fillthis = "";
2193 
2194  //we loop over all tracks in the event
2195  for (Int_t i=0; i<aodEvent->GetNumberOfTracks(); i++){
2196  AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(i));
2197  if(!aodTrack) AliFatal("Not a standard AOD");
2198 
2199  //quick quality cut
2200  if(aodTrack->GetITSNcls() < 1) continue;
2201  if(aodTrack->GetTPCNcls() < 1) continue;
2202  if (aodTrack->GetStatus()&AliESDtrack::kITSpureSA) continue;
2203  if (!(aodTrack->GetStatus()&AliESDtrack::kITSin)) continue;
2204  if (aodTrack->GetID() < 0) continue;
2205  Double_t covtest[21];
2206  if (!aodTrack->GetCovarianceXYZPxPyPz(covtest)) continue;
2207 
2208 
2209  Int_t mcLabelParticle = -1;
2210  Int_t pdgParticle = -1;
2211  mcLabelParticle = aodTrack->GetLabel();
2212 
2213  numberofparticles++;
2214 
2215  // we fill histograms with information of the track
2216  Double_t pt_track = aodTrack->Pt();
2217  Double_t momentum_track = aodTrack->P();
2218  Int_t numberOfITS = aodTrack->GetITSNcls();
2219  Int_t numberOfTPC = aodTrack->GetTPCNcls();
2220  Double_t nSigmaTPC = 0;
2221  Double_t nSigmaTOF = 0;
2222  Int_t pionPIDnumber = 2;
2223  Int_t kaonPIDnumber = 3;
2224  Int_t TPCok = 0;
2225  Int_t TOFok = 0;
2226 
2227  AliAODPidHF* trackPIDHF = (AliAODPidHF*)fCuts->GetPidHF();
2228  TPCok = trackPIDHF->GetnSigmaTPC(aodTrack, pionPIDnumber, nSigmaTPC);
2229  TOFok = trackPIDHF->GetnSigmaTOF(aodTrack, pionPIDnumber, nSigmaTOF);
2230 
2231  AliExternalTrackParam particleTrack;
2232  particleTrack.CopyFromVTrack(aodTrack);
2233  Double_t d0[2],covd0[3];
2234  particleTrack.PropagateToDCA(primaryVertex,bz,100.,d0,covd0);
2235 
2236 
2237  //we check if the particle is a signal track
2238  Bool_t isDesiredCandidate = kFALSE;
2239  if(fUseMCInfo){
2240  TMatrix &particleMatrix = *BPlustoD0PiLabelMatrix;
2241  for (Int_t k = 0; k < BPlustoD0PiLabelMatrix->GetNrows(); ++k){
2242  if(mcLabelParticle == (Int_t)particleMatrix(k,0)){
2243  isDesiredCandidate = kTRUE;
2244  break;
2245  }
2246  }
2247  }
2248 
2249  if(fUseMCInfo){
2250  if(IsTrackInjected(aodTrack,header,mcTrackArray) && !isDesiredCandidate && fQuickSignalAnalysis == 2) continue;
2251  }
2252 
2253  Int_t daughterType = 2;
2254 
2255  Int_t histType = 0;
2256  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
2257  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
2258  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
2259  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
2260 
2261  for (Int_t j = 0; j < 10; ++j)
2262  {
2263  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
2264 
2265  }
2266 
2267  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
2268  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
2269  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
2270  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
2271 
2272  if(isDesiredCandidate)
2273  {
2274  histType = 1;
2275  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
2276  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
2277  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
2278  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
2279 
2280  for (Int_t j = 0; j < 10; ++j)
2281  {
2282  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
2283 
2284  }
2285 
2286  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
2287  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
2288  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
2289  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
2290  }
2291 
2292 
2293  //we apply a number of cuts on the particle
2294  Bool_t bCut = kFALSE;
2295 
2296 
2297  //we apply a PID cut, 2 is used to indicate we look for a pion
2298  if(!(fCuts->SelectPID(aodTrack,2))){
2299  if(isDesiredCandidate) {
2300  ((TH1F*)fDaughterHistogramArrayExtra[2][1])->Fill(2);
2301  } else ((TH1F*)fDaughterHistogramArrayExtra[2][0])->Fill(2);
2302  bCut = kTRUE;
2303  }
2304 
2305  if(aodTrack->GetITSNcls() < fCuts->GetMinITSNclsBPlusPion()){
2306  if(isDesiredCandidate) {
2307  ((TH1F*)fDaughterHistogramArrayExtra[2][1])->Fill(3);
2308  } else ((TH1F*)fDaughterHistogramArrayExtra[2][0])->Fill(3);
2309  bCut = kTRUE;
2310  }
2311 
2312  if(aodTrack->GetTPCNcls() < fCuts->GetMinTPCNclsBPlusPion()){
2313  if(isDesiredCandidate) {
2314  ((TH1F*)fDaughterHistogramArrayExtra[2][1])->Fill(4);
2315  } else ((TH1F*)fDaughterHistogramArrayExtra[2][0])->Fill(4);
2316  bCut = kTRUE;
2317  }
2318 
2319  if(fCuts->UseITSRefitBPlusPion()==kTRUE){
2320  if(!(aodTrack->GetStatus()&AliESDtrack::kITSrefit)) {
2321  if(isDesiredCandidate) {
2322  ((TH1F*)fDaughterHistogramArrayExtra[2][1])->Fill(5);
2323  } else ((TH1F*)fDaughterHistogramArrayExtra[2][0])->Fill(5);
2324  bCut = kTRUE;
2325  }
2326  }
2327 
2328  if(fCuts->UseTPCRefitBPlusPion()==kTRUE){
2329  if((!(aodTrack->GetStatus()&AliESDtrack::kTPCrefit))) {
2330  if(isDesiredCandidate) {
2331  ((TH1F*)fDaughterHistogramArrayExtra[2][1])->Fill(6);
2332  } else ((TH1F*)fDaughterHistogramArrayExtra[2][0])->Fill(6);
2333  bCut = kTRUE;
2334  }
2335  }
2336 
2337  if(fCuts->UseFilterBitBPlusPion()==kTRUE){
2338  if(!(aodTrack->TestFilterMask(BIT(fCuts->GetFilterBitBPlusPion())))) {
2339  if(isDesiredCandidate) {
2340  ((TH1F*)fDaughterHistogramArrayExtra[2][1])->Fill(7);
2341  } else ((TH1F*)fDaughterHistogramArrayExtra[2][0])->Fill(7);
2342  bCut = kTRUE;
2343  }
2344  }
2345 
2346 
2347  if(aodTrack->Pt() < fCuts->GetMinPtBPlusPion()){
2348  if(isDesiredCandidate) {
2349  ((TH1F*)fDaughterHistogramArrayExtra[2][1])->Fill(8);
2350  } else ((TH1F*)fDaughterHistogramArrayExtra[2][0])->Fill(8);
2351  bCut = kTRUE;
2352  }
2353 
2354 
2355  if(TMath::Abs(d0[0]) < fCuts->GetMind0BPlusPion()){
2356  if(isDesiredCandidate) {
2357  ((TH1F*)fDaughterHistogramArrayExtra[2][1])->Fill(12);
2358  } else ((TH1F*)fDaughterHistogramArrayExtra[2][0])->Fill(12);
2359  bCut = kTRUE;
2360  }
2361 
2362  if(TMath::Abs(aodTrack->Eta()) > fCuts->GetMaxAbsEtaBPlusPion()){
2363  if(isDesiredCandidate) {
2364  ((TH1F*)fDaughterHistogramArrayExtra[2][1])->Fill(9);
2365  } else ((TH1F*)fDaughterHistogramArrayExtra[2][0])->Fill(9);
2366  bCut = kTRUE;
2367  }
2368 
2369  Bool_t bHardSelectionArrayITS[7] = {kFALSE};
2370  fCuts->GetHardSelectionArrayITSBPlusPion(bHardSelectionArrayITS);
2371  Bool_t bSoftSelectionArrayITS[7] = {kFALSE};
2372  fCuts->GetSoftSelectionArrayITSBPlusPion(bSoftSelectionArrayITS);
2373 
2374  Bool_t bHardITSPass = kTRUE;
2375  for (Int_t j = 0; j < 7; ++j)
2376  {
2377  if(bHardSelectionArrayITS[j])
2378  {
2379  if(!aodTrack->HasPointOnITSLayer(j)) bHardITSPass = kFALSE;
2380  }
2381  }
2382 
2383  Int_t nCounterSoftSelection = 0;
2384  Bool_t bSoftITSPass = kTRUE;
2385  for (Int_t j = 0; j < 7; ++j)
2386  {
2387  if(bSoftSelectionArrayITS[j])
2388  {
2389  if(aodTrack->HasPointOnITSLayer(j)) nCounterSoftSelection++;
2390  }
2391  }
2392  if(nCounterSoftSelection < fCuts->GetNSoftITSCutBPlusPion()) bSoftITSPass = kFALSE;
2393 
2394  if(!bHardITSPass){
2395  if(isDesiredCandidate) {
2396  ((TH1F*)fDaughterHistogramArrayExtra[2][1])->Fill(10);
2397  } else ((TH1F*)fDaughterHistogramArrayExtra[2][0])->Fill(10);
2398  bCut = kTRUE;
2399  }
2400 
2401  if(!bSoftITSPass){
2402  if(isDesiredCandidate) {
2403  ((TH1F*)fDaughterHistogramArrayExtra[2][1])->Fill(11);
2404  } else ((TH1F*)fDaughterHistogramArrayExtra[2][0])->Fill(11);
2405  bCut = kTRUE;
2406  }
2407 
2408 
2409  if(!isDesiredCandidate && fQuickSignalAnalysis == 1) bCut = kTRUE;
2410 
2411  if(bCut) {
2412  if(isDesiredCandidate) {
2413  ((TH1F*)fDaughterHistogramArrayExtra[2][1])->Fill(0);
2414  } else ((TH1F*)fDaughterHistogramArrayExtra[2][0])->Fill(0);
2415  continue;
2416  }
2417 
2418  //we fill histograms with track information of the tracks that pass the cuts
2419  histType = 2;
2420  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
2421  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
2422  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
2423  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
2424 
2425  for (Int_t j = 0; j < 10; ++j)
2426  {
2427  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
2428 
2429  }
2430 
2431  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
2432  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
2433  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
2434  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
2435 
2436  if(isDesiredCandidate)
2437  {
2438  histType = 3;
2439  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
2440  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
2441  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
2442  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
2443 
2444  for (Int_t j = 0; j < 10; ++j)
2445  {
2446  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
2447 
2448  }
2449 
2450  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
2451  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
2452  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
2453  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
2454  }
2455 
2456  fBPlusPionTracks->push_back(i);
2457  numberofparticlesused++;
2458  }
2459 
2460  // cout << "BPlus pions used: " << numberofparticlesused << endl;
2461 
2462  ((TH1F*)fDaughterHistogramArray[2][0][10])->Fill(numberofparticles);
2463  ((TH1F*)fDaughterHistogramArray[2][1][10])->Fill(numberofparticlesused);
2464  return;
2465 }
2466 //-------------------------------------------------------------------------------------
2467 void AliAnalysisTaskSEBPlustoD0Pi::D0Selection(AliAODEvent* aodEvent, AliAODVertex *primaryVertex, Double_t bz,TClonesArray * mcTrackArray, TMatrix * BPlustoD0PiLabelMatrix, TClonesArray * D0TracksFromFriendFile, AliAODMCHeader * header){
2468 
2469  TString fillthis = "";
2470 
2472 
2473  //next we loop over all the D0 candidates
2474  for (Int_t j = 0; j < D0TracksFromFriendFile->GetEntriesFast(); j++)
2475  {
2476 
2477  //we get the track of the D0
2478  AliAODRecoDecayHF2Prong * trackD0 = (AliAODRecoDecayHF2Prong*)(D0TracksFromFriendFile->At(j));
2479  if(!trackD0) {std::cout << "found none" << std::endl; continue;}
2480  if(trackD0 == nullptr) {std::cout << "found nullptr" << std::endl; continue;}
2481 
2482  if(!(vHF->FillRecoCand(aodEvent,trackD0))) //Fill the data members of the candidate only if they are empty.
2483  {
2484  fCEvents->Fill(12); //monitor how often this fails
2485  continue;
2486  }
2487 
2488  AliAODTrack * trackFirstDaughter = (AliAODTrack*)(trackD0->GetDaughter(0));
2489  AliAODTrack * trackSecondDaughter = (AliAODTrack*)(trackD0->GetDaughter(1));
2490  if(!D0FirstDaughterSelection(trackFirstDaughter, primaryVertex, bz, mcTrackArray, BPlustoD0PiLabelMatrix,header)) continue;
2491  if(!D0SecondDaughterSelection(trackSecondDaughter, primaryVertex, bz, mcTrackArray, BPlustoD0PiLabelMatrix,header)) continue;
2492 
2493 
2494  AliAODVertex *vertexMother = (AliAODVertex*)trackD0->GetSecondaryVtx();
2495 
2496  //we save the pdgcode of the used particle and its mother and check if it is a desired candidate
2497  Int_t pdgCodeMother = -1;
2498  Float_t pdgCodeGrandMother = -1;
2499  Bool_t isDesiredCandidate = kFALSE;
2500  Int_t motherType, histType;
2501  motherType = 0;
2502  Int_t mcLabelD0 = -1;
2503 
2504  if(fUseMCInfo)
2505  {
2506  mcLabelD0 = MatchCandidateToMonteCarlo(421,trackD0,mcTrackArray,BPlustoD0PiLabelMatrix);
2507 
2508  if(mcLabelD0 >= 0)
2509  {
2510  isDesiredCandidate = kTRUE;
2511 
2512  Int_t mcLabelFirstTrack = -1;
2513  mcLabelFirstTrack = trackFirstDaughter->GetLabel();
2514 
2515  if(mcLabelFirstTrack >= 0)
2516  {
2517  AliAODMCParticle *mcParticleFirstTrack = (AliAODMCParticle*)mcTrackArray->At(mcLabelFirstTrack);
2518  AliAODMCParticle *mcMotherParticle = (AliAODMCParticle*)mcTrackArray->At(mcLabelD0);
2519 
2520  if(mcParticleFirstTrack && mcMotherParticle)
2521  {
2522  pdgCodeMother = mcMotherParticle->GetPdgCode();
2523 
2524  Double_t vertex_distance = TMath::Sqrt((vertexMother->GetX() - mcParticleFirstTrack->Xv())*(vertexMother->GetX() - mcParticleFirstTrack->Xv()) + (vertexMother->GetY() - mcParticleFirstTrack->Yv())*(vertexMother->GetY() - mcParticleFirstTrack->Yv()) + (vertexMother->GetZ() - mcParticleFirstTrack->Zv())*(vertexMother->GetZ() - mcParticleFirstTrack->Zv()));
2525  ((TH1F*)fMotherHistogramArrayExtra[motherType][4])->Fill(vertex_distance);
2526 
2527  Double_t momentum_resolution = TMath::Sqrt((trackD0->Px() - mcMotherParticle->Px())*(trackD0->Px() - mcMotherParticle->Px()) + (trackD0->Py() - mcMotherParticle->Py())*(trackD0->Py() - mcMotherParticle->Py()) + (trackD0->Pz() - mcMotherParticle->Pz())*(trackD0->Pz() - mcMotherParticle->Pz()));
2528  ((TH1F*)fMotherHistogramArrayExtra[motherType][6])->Fill(momentum_resolution);
2529  }
2530  }
2531  }
2532  }
2533 
2534  // We fill the histograms
2535  histType = 0;
2536  FillD0Histograms(trackD0, primaryVertex, bz, motherType, histType);
2537  if(isDesiredCandidate && fUseMCInfo){
2538  histType = 1;
2539  FillD0Histograms(trackD0, primaryVertex, bz, motherType, histType,pdgCodeMother);
2540  }
2541 
2542  // Here we apply cuts on the particle
2543  Bool_t cutMother = kFALSE;
2544 
2545  Bool_t bCutArray[29] = {0};
2546  Int_t cutReturnValue = fCuts->IsD0forD0ptbinSelected(trackD0, 0, aodEvent, bCutArray);
2547  if(cutReturnValue == -1) cutMother = kTRUE;
2548  if(cutReturnValue == 0) cutMother = kTRUE;
2549 
2550 
2551  if(fGetCutInfo == kTRUE)
2552  {
2553  for (Int_t k = 0; k < 29; ++k)
2554  {
2555  if (bCutArray[k] == kTRUE){
2556  if(isDesiredCandidate){
2557  ((TH1F*)fMotherHistogramArrayExtra[motherType][1])->Fill(k+1);
2558  } else ((TH1F*)fMotherHistogramArrayExtra[motherType][0])->Fill(k+1);
2559  cutMother = kTRUE;
2560  }
2561  }
2562  }
2563 
2564 
2565  if(!isDesiredCandidate && fQuickSignalAnalysis == 1) cutMother = kTRUE;
2566 
2567  if(cutMother){
2568  if(isDesiredCandidate){
2569  ((TH1F*)fMotherHistogramArrayExtra[motherType][1])->Fill(0);
2570  } else ((TH1F*)fMotherHistogramArrayExtra[motherType][0])->Fill(0);
2571  // delete vertexMother; vertexMother = nullptr;
2572  // delete trackD0; trackD0 = nullptr;
2573  continue;
2574  }
2575 
2576  // We fill the cut histograms
2577  histType = 2;
2578  FillD0Histograms(trackD0, primaryVertex, bz, motherType, histType);
2579  if(isDesiredCandidate && fUseMCInfo){
2580  histType = 3;
2581  FillD0Histograms(trackD0, primaryVertex, bz, motherType, histType,pdgCodeMother);
2582  }
2583 
2584  //we save the location of the D0 candidate
2585  fD0Tracks->push_back(j);
2586  }
2587 
2588  delete vHF; vHF = nullptr;
2589  return;
2590 }
2591 //-------------------------------------------------------------------------------------
2592 void AliAnalysisTaskSEBPlustoD0Pi::BPlusSelection(AliAODEvent* aodEvent, AliAODVertex *primaryVertex, Double_t bz, TClonesArray * mcTrackArray, TMatrix * BPlustoD0PiLabelMatrix, TClonesArray * D0TracksFromFriendFile, AliAODMCHeader * header) {
2593 
2594  // Int_t iClonesArray = 0;
2595 
2596  TString fillthis = "";
2597 
2598  //we loop over all the D0 candidates
2599  for (Int_t j = 0; j < (Int_t)fD0Tracks->size(); j++)
2600  {
2601 
2602  //Save current Object count
2603  Int_t ObjectNumber = TProcessID::GetObjectCount();
2604 
2605  //we get the track of the D0
2606  AliAODRecoDecayHF2Prong * trackD0 = (AliAODRecoDecayHF2Prong*)(D0TracksFromFriendFile->At(fD0Tracks->at(j)));
2607  if(!trackD0) {std::cout << "found none" << std::endl; continue;}
2608  if(trackD0 == nullptr) {std::cout << "found nullptr" << std::endl; continue;}
2609 
2610  //we loop over all the BPlus pion candidates
2611  for (Int_t i = 0; i < (Int_t)fBPlusPionTracks->size(); i++)
2612  {
2613 
2614  //we get the track of the BPlus daughter
2615  AliAODTrack * trackBPlusPion = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(fBPlusPionTracks->at(i)));
2616  if (!trackBPlusPion) continue;
2617 
2618  //we check if the IDs of the tracks are different
2619  AliAODTrack* twoProngdaughter0 = (AliAODTrack*)trackD0->GetDaughter(0);
2620  AliAODTrack* twoProngdaughter1 = (AliAODTrack*)trackD0->GetDaughter(1);
2621  UShort_t idProng0 = twoProngdaughter0->GetID();
2622  UShort_t idProng1 = twoProngdaughter1->GetID();
2623 
2624  if (trackBPlusPion->GetID() == idProng0 || trackBPlusPion->GetID() == idProng1) continue;
2625 
2626  UInt_t prongsD0[2];
2627  prongsD0[0] = 211;
2628  prongsD0[1] = 321;
2629 
2630 
2631  UInt_t prongsD02[2];
2632  prongsD02[1] = 211;
2633  prongsD02[0] = 321;
2634 
2635  // D0 window - invariant mass cut
2636  Double_t invariantMassD0 = trackD0->InvMass(2, prongsD0);
2637  Double_t invariantMassD02 = trackD0->InvMass(2, prongsD02);
2638 
2639  Double_t pdgMassD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2640  Double_t massWindowD0 = 0.06; //GeV/c^2 //-----------------------------------------fcuts get mass window --------------------------------------------------------------------------
2641  Int_t pdgD0 = 421;
2642 
2643  //we check if the pions have the opposite charge
2644  Bool_t bWrongSign = kFALSE;
2645  if (trackBPlusPion->Charge() == -1)
2646  {
2647  if ((fCuts->SelectPID(((AliAODTrack*)trackD0->GetDaughter(0)), 2)) && (fCuts->SelectPID(((AliAODTrack*)trackD0->GetDaughter(1)), 3))) bWrongSign = kFALSE;
2648  else if ((fCuts->SelectPID(((AliAODTrack*)trackD0->GetDaughter(0)), 3)) && (fCuts->SelectPID(((AliAODTrack*)trackD0->GetDaughter(1)), 2))) bWrongSign = kTRUE;
2649  else continue;
2650  if (TMath::Abs(invariantMassD0 - pdgMassD0) > massWindowD0) continue;
2651 
2652  } else if (trackBPlusPion->Charge() == 1) {
2653  pdgD0 = -421;
2654  if ((fCuts->SelectPID(((AliAODTrack*)trackD0->GetDaughter(0)), 3)) && (fCuts->SelectPID(((AliAODTrack*)trackD0->GetDaughter(1)), 2))) bWrongSign = kFALSE;
2655  else if ((fCuts->SelectPID(((AliAODTrack*)trackD0->GetDaughter(0)), 2)) && (fCuts->SelectPID(((AliAODTrack*)trackD0->GetDaughter(1)), 3))) bWrongSign = kTRUE;
2656  else continue;
2657  if (TMath::Abs(invariantMassD02 - pdgMassD0) > massWindowD0) continue;
2658  }
2659 
2660  //location BPlus pion rotation around PV
2661  for (Int_t iRot = 0; iRot < fNumberOfRotations + 1; ++iRot)
2662  {
2663  //we create a copy of the track that we will rotate
2664  AliAODTrack * trackBPlusPionRotated = new AliAODTrack(*trackBPlusPion);
2665 
2666  //for iRot == 0, we use the original unrotated track. For iRot > 0 we rotate the track and set the label to -1
2667  if (iRot != 0)
2668  {
2669  //should still check if track is already at PV
2670  Double_t dPhiRotated = trackBPlusPionRotated->Phi() + TMath::Pi() - (TMath::Pi() * fDegreePerRotation * fNumberOfRotations / (180.0 * 2.0)) + (TMath::Pi() * fDegreePerRotation * iRot / 180.0);
2671  trackBPlusPionRotated->SetPhi(dPhiRotated);
2672  }
2673 
2674 
2675  //we use the BPlus pion and D0 tracks to reconstruct the vertex for the BPlus
2676  AliExternalTrackParam firstTrack;
2677  firstTrack.CopyFromVTrack(trackBPlusPionRotated);
2678  AliExternalTrackParam secondTrack;
2679  secondTrack.CopyFromVTrack(trackD0);
2680 
2681  UInt_t prongs[2];
2682  prongs[1] = 421;
2683  prongs[0] = 211;
2684 
2685  // we calculate the vertex of the mother candidate
2686  TObjArray daughterTracks;
2687 
2688  daughterTracks.Add(&firstTrack);
2689  daughterTracks.Add(&secondTrack);
2690  Double_t dispersion = 0;
2691  AliAODVertex *vertexMother = RecalculateVertex(primaryVertex, &daughterTracks, bz, dispersion);
2692  if (!vertexMother)
2693  {
2694  delete vertexMother; vertexMother = nullptr;
2695  delete trackBPlusPionRotated; trackBPlusPionRotated = nullptr;
2696  continue;
2697  }
2698 
2699 
2700  //use the new vertex to create the BPlus candidate
2701  Double_t xdummy = 0., ydummy = 0., dca, e[2];
2702  Double_t d0z0[2], covd0z0[3], d0[2], d0err[2];
2703 
2704 
2705  firstTrack.PropagateToDCA(vertexMother, bz, 100., d0z0, covd0z0);
2706  secondTrack.PropagateToDCA(vertexMother, bz, 100., d0z0, covd0z0);
2707 
2708  //we reconstruct the mother decay prong
2709  Double_t px[2], py[2], pz[2];
2710  px[0] = firstTrack.Px();
2711  py[0] = firstTrack.Py();
2712  pz[0] = firstTrack.Pz();
2713  px[1] = secondTrack.Px();
2714  py[1] = secondTrack.Py();
2715  pz[1] = secondTrack.Pz();
2716 
2717  UShort_t id[2];
2718  id[0] = firstTrack.GetID();
2719  id[1] = secondTrack.GetID();
2720 
2721  firstTrack.PropagateToDCA(primaryVertex, bz, 100., d0z0, covd0z0);
2722  d0[0] = d0z0[0];
2723  d0err[0] = TMath::Sqrt(covd0z0[0]);
2724  secondTrack.PropagateToDCA(primaryVertex, bz, 100., d0z0, covd0z0);
2725  d0[1] = d0z0[0];
2726  d0err[1] = TMath::Sqrt(covd0z0[0]);
2727 
2728 
2729  dca = secondTrack.GetDCA(&firstTrack, bz, xdummy, ydummy);
2730 
2731 
2732  Short_t chargeMother = trackD0->Charge() + trackBPlusPionRotated->Charge();
2733  // cout << " charge: " << chargeMother << endl;
2734  AliAODRecoDecayHF2Prong trackBPlus(vertexMother, px, py, pz, d0, d0err, dca);
2735 
2736  trackBPlus.SetCharge(chargeMother);
2737  trackBPlus.GetSecondaryVtx()->AddDaughter(trackBPlusPionRotated);
2738  trackBPlus.GetSecondaryVtx()->AddDaughter(trackD0);
2739  trackBPlus.SetPrimaryVtxRef((AliAODVertex*)aodEvent->GetPrimaryVertex());
2740  trackBPlus.SetProngIDs(2, id);
2741 
2742  // Fiducial cut
2743  if(TMath::Abs(trackBPlus.Y(521)) > 0.8) {
2744  delete vertexMother; vertexMother = nullptr;
2745  delete trackBPlusPionRotated; trackBPlusPionRotated = nullptr;
2746  continue;
2747  }
2748 
2749  // We check if the signal is injected, optionally we can reject injected signals
2750  Bool_t fCheckInjected = kTRUE; //temp
2751  Bool_t fRemoveInjected = kFALSE; //temp
2752  Bool_t bIsInjected = kFALSE;
2753  if(fUseMCInfo){
2754  if(fCheckInjected) bIsInjected = IsCandidateInjected(&trackBPlus, header,mcTrackArray);
2755  if(fCheckInjected && fRemoveInjected && bIsInjected) {
2756  delete vertexMother; vertexMother = nullptr;
2757  delete trackBPlusPionRotated; trackBPlusPionRotated = nullptr;
2758  continue;
2759  }
2760  }
2761 
2762  // We check if the BPlus candidate is a true signal in Monte Carlo
2763  Bool_t isDesiredCandidate = kFALSE;
2764  Int_t mcLabelBPlus = -1;
2765  fillthis = "";
2766  Int_t motherType, histType;
2767  motherType = 1;
2768 
2769  if(fUseMCInfo)
2770  {
2771  mcLabelBPlus = MatchCandidateToMonteCarlo(521,&trackBPlus,mcTrackArray,BPlustoD0PiLabelMatrix);
2772  if (mcLabelBPlus >= 0 && trackBPlusPionRotated->GetLabel() >= 0 && iRot == 0) isDesiredCandidate = kTRUE;
2773  }
2774 
2775  Double_t ptMother = trackBPlus.Pt();
2776 
2777  histType = 0;
2778 
2779  if (isDesiredCandidate)
2780  {
2781  AliAODMCParticle *mcTrackBPlusPion = (AliAODMCParticle*)mcTrackArray->At(trackBPlusPion->GetLabel());
2782  AliAODMCParticle *mcTrackBPlus = (AliAODMCParticle*)mcTrackArray->At(mcLabelBPlus);
2783 
2784  Double_t vertex_distance = TMath::Sqrt((primaryVertex->GetX() - mcTrackBPlusPion->Xv()) * (primaryVertex->GetX() - mcTrackBPlus->Xv()) + (primaryVertex->GetY() - mcTrackBPlus->Yv()) * (primaryVertex->GetY() - mcTrackBPlus->Yv()) + (primaryVertex->GetZ() - mcTrackBPlus->Zv()) * (primaryVertex->GetZ() - mcTrackBPlus->Zv()));
2785  ((TH1F*)fMotherHistogramArrayExtra[motherType][4])->Fill(vertex_distance);
2786 
2787  Double_t momentum_resolution = TMath::Sqrt((trackBPlus.Px() - mcTrackBPlus->Px()) * (trackBPlus.Px() - mcTrackBPlus->Px()) + (trackBPlus.Py() - mcTrackBPlus->Py()) * (trackBPlus.Py() - mcTrackBPlus->Py()) + (trackBPlus.Pz() - mcTrackBPlus->Pz()) * (trackBPlus.Pz() - mcTrackBPlus->Pz()));
2788  ((TH1F*)fMotherHistogramArrayExtra[motherType][6])->Fill(momentum_resolution);
2789 
2790  }
2791 
2792 
2793  if (!bWrongSign)
2794  {
2795  // We fill the histograms
2796  histType = 0;
2797  FillBPlusHistograms(&trackBPlus, primaryVertex, bz, motherType, histType);
2798 
2799  if (isDesiredCandidate)
2800  {
2801  histType = 1;
2802  FillBPlusHistograms(&trackBPlus, primaryVertex, bz, motherType, histType);
2803  }
2804  }
2805 
2806 
2807  // We apply cuts
2808  Bool_t cutMother = kFALSE;
2809 
2810  Bool_t bCutArray[68] = {0};
2811  Int_t numberOfCuts = 68;
2812  Int_t cutReturnValue = fCuts->IsSelected(&trackBPlus, 0, aodEvent, bCutArray);
2813  if(cutReturnValue == -1) cutMother = kTRUE;
2814  if(cutReturnValue == 0) cutMother = kTRUE;
2815 
2816 
2817  // We save information about the cuts
2818  TString histName = "";
2819  Double_t invariantMassMother = trackBPlus.InvMass(2,prongs);
2820  Double_t pdgMassMother=TDatabasePDG::Instance()->GetParticle(521)->Mass();
2821  Double_t massWindow = fHistMassWindow; //GeV/c^2
2822  if(fGetCutInfo == kTRUE)
2823  {
2824  for (Int_t n = 0; n < 68; ++n)
2825  {
2826  if(bCutArray[n] == kTRUE){
2827  if(isDesiredCandidate){
2828  ((TH1F*)fMotherHistogramArrayExtra[motherType][1])->Fill(n+1);
2829  } else ((TH1F*)fMotherHistogramArrayExtra[motherType][0])->Fill(n+1);
2830  cutMother = kTRUE;
2831  }
2832  }
2833 
2834  if (TMath::Abs(invariantMassMother-pdgMassMother)<massWindow){
2835  for (Int_t l = 0; l < numberOfCuts; ++l) //total
2836  {
2837  if(bCutArray[l] == kFALSE) continue;
2838  for (Int_t j = 0; j < numberOfCuts; ++j)
2839  {
2840  if(bCutArray[j] == kFALSE) continue;
2841  if(isDesiredCandidate == kFALSE) histName ="cutEffectBackground";
2842  if(isDesiredCandidate == kTRUE) histName ="cutEffectSignal";
2843  ((TH2I*)(fOutputBPlusMC->FindObject(histName)))->Fill(l,j);
2844  }
2845  }
2846 
2847  for (Int_t l = 0; l < numberOfCuts; ++l) //unique
2848  {
2849  if(bCutArray[l] == kFALSE) continue;
2850  Bool_t bFill = kTRUE;
2851  for (Int_t j = 0; j < numberOfCuts; ++j)
2852  {
2853  if(l==j) continue;
2854  if(bCutArray[j] == kTRUE)
2855  {
2856  bFill = kFALSE;
2857  break;
2858  }
2859 
2860  }
2861  if(bFill == kTRUE)
2862  {
2863  if(isDesiredCandidate == kFALSE) histName ="cutEffectUniqueBackground";
2864  if(isDesiredCandidate == kTRUE) histName ="cutEffectUniqueSignal";
2865  ((TH1I*)(fOutputBPlusMC->FindObject(histName)))->Fill(l);
2866  }
2867  }
2868  }
2869  }
2870 
2871 
2872  if(!isDesiredCandidate && fQuickSignalAnalysis == 1) cutMother = kTRUE;
2873 
2874  if(cutMother)
2875  {
2876  if(isDesiredCandidate)
2877  {
2878  ((TH1F*)fMotherHistogramArrayExtra[motherType][1])->Fill(0);
2879  } else ((TH1F*)fMotherHistogramArrayExtra[motherType][0])->Fill(0);
2880  delete vertexMother; vertexMother = nullptr;
2881  delete trackBPlusPionRotated; trackBPlusPionRotated = nullptr;
2882  continue;
2883  }
2884 
2885  // isDesiredCandidate = kFALSE;
2886 
2887  // if (4.9 < invariantMassMother && invariantMassMother < 5.2)
2888  // {
2889  // if(!bWrongSign)
2890  // {
2891  // Int_t mcLabelBPlusPion = trackBPlusPion->GetLabel();
2892  // Int_t mcLabelD0first = ((AliAODTrack*)trackD0->GetDaughter(0))->GetLabel();
2893  // Int_t mcLabelD0second = ((AliAODTrack*)trackD0->GetDaughter(1))->GetLabel();
2894 
2895  // AliAODMCParticle * mcBPlusPion = (AliAODMCParticle*)mcTrackArray->At(mcLabelBPlusPion);
2896  // AliAODMCParticle * mcD0first = (AliAODMCParticle*)mcTrackArray->At(mcLabelD0first);
2897  // AliAODMCParticle * mcD0second = (AliAODMCParticle*)mcTrackArray->At(mcLabelD0second);
2898 
2899  // if(mcBPlusPion)
2900  // {
2901  // while(mcBPlusPion->GetMother() >= 0)
2902  // {
2903  // mcBPlusPion = (AliAODMCParticle*)mcTrackArray->At(mcBPlusPion->GetMother());
2904  // fillthis="particle_pdgBPlusPion";
2905  // ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(TMath::Abs(mcBPlusPion->GetPdgCode()));
2906  // }
2907  // }
2908 
2909  // if(mcD0first)
2910  // {
2911  // while(mcD0first->GetMother() >= 0)
2912  // {
2913  // mcD0first = (AliAODMCParticle*)mcTrackArray->At(mcD0first->GetMother());
2914  // fillthis="particle_pdgD0First";
2915  // ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(TMath::Abs(mcD0first->GetPdgCode()));
2916  // }
2917  // }
2918 
2919  // if(mcD0second)
2920  // {
2921  // while(mcD0second->GetMother() >= 0)
2922  // {
2923  // mcD0second = (AliAODMCParticle*)mcTrackArray->At(mcD0second->GetMother());
2924  // fillthis="particle_pdgD0Second";
2925  // ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(TMath::Abs(mcD0second->GetPdgCode()));
2926  // }
2927  // }
2928 
2929  // mcBPlusPion = (AliAODMCParticle*)mcTrackArray->At(mcLabelBPlusPion);
2930  // mcD0first = (AliAODMCParticle*)mcTrackArray->At(mcLabelD0first);
2931  // mcD0second = (AliAODMCParticle*)mcTrackArray->At(mcLabelD0second);
2932 
2933  // if(mcBPlusPion && mcD0first && mcD0second)
2934  // {
2935  // if(mcD0first->GetMother() == mcD0second->GetMother() && mcD0first->GetMother() >= 0)
2936  // {
2937  // AliAODMCParticle * D0Mother = (AliAODMCParticle*)mcTrackArray->At(mcD0first->GetMother());
2938  // if(D0Mother->GetMother() == mcBPlusPion->GetMother() && D0Mother->GetMother() >= 0)
2939  // {
2940  // AliAODMCParticle * finalMother = (AliAODMCParticle*)mcTrackArray->At(D0Mother->GetMother());
2941  // fillthis="particle_pdgAll";
2942  // ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(TMath::Abs(finalMother->GetPdgCode()));
2943  // fillthis="particle_pdgAllInvMass";
2944  // ((TH2F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(TMath::Abs(finalMother->GetPdgCode()),invariantMassMother);
2945  // if(finalMother->GetPdgCode()==511)
2946  // {
2947  // // isDesiredCandidate = kTRUE;
2948  // for (Int_t iDaughter = 0; iDaughter < finalMother->GetNDaughters(); ++iDaughter) //will work up to ten daughters
2949  // {
2950  // AliAODMCParticle* daughterfinalMother = (AliAODMCParticle*)mcTrackArray->At(finalMother->GetDaughter(0)+iDaughter);
2951  // fillthis="particle_daughterPdgOneStep511";
2952  // ((TH2F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(iDaughter,TMath::Abs(daughterfinalMother->GetPdgCode()));
2953  // }
2954  // for (Int_t iDaughter = 0; iDaughter < D0Mother->GetNDaughters(); ++iDaughter) //will work up to ten daughters
2955  // {
2956  // AliAODMCParticle* daughterD0Mother = (AliAODMCParticle*)mcTrackArray->At(D0Mother->GetDaughter(0)+iDaughter);
2957  // fillthis="particle_daughterPdgOneStep511";
2958  // ((TH2F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(iDaughter+10,TMath::Abs(daughterD0Mother->GetPdgCode()));
2959  // }
2960  // }
2961  // if(finalMother->GetPdgCode()==521)
2962  // {
2963  // // isDesiredCandidate = kTRUE;
2964  // for (Int_t iDaughter = 0; iDaughter < finalMother->GetNDaughters(); ++iDaughter) //will work up to ten daughters
2965  // {
2966  // AliAODMCParticle* daughterfinalMother = (AliAODMCParticle*)mcTrackArray->At(finalMother->GetDaughter(0)+iDaughter);
2967  // fillthis="particle_daughterPdgOneStep521";
2968  // ((TH2F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(iDaughter,TMath::Abs(daughterfinalMother->GetPdgCode()));
2969  // }
2970  // for (Int_t iDaughter = 0; iDaughter < D0Mother->GetNDaughters(); ++iDaughter) //will work up to ten daughters
2971  // {
2972  // AliAODMCParticle* daughterD0Mother = (AliAODMCParticle*)mcTrackArray->At(D0Mother->GetDaughter(0)+iDaughter);
2973  // fillthis="particle_daughterPdgOneStep521";
2974  // ((TH2F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(iDaughter+10,TMath::Abs(daughterD0Mother->GetPdgCode()));
2975  // }
2976  // }
2977  // }
2978  // }
2979 
2980  // if(mcD0first->GetMother() == mcD0second->GetMother() && mcD0first->GetMother() >= 0)
2981  // {
2982  // AliAODMCParticle * D0Mother = (AliAODMCParticle*)mcTrackArray->At(mcD0first->GetMother());
2983  // AliAODMCParticle * D0GrandMother = (AliAODMCParticle*)mcTrackArray->At(D0Mother->GetMother());
2984  // if(D0GrandMother->GetMother() == mcBPlusPion->GetMother() && D0GrandMother->GetMother() >= 0)
2985  // {
2986  // AliAODMCParticle * finalMother = (AliAODMCParticle*)mcTrackArray->At(D0GrandMother->GetMother());
2987  // fillthis="particle_pdgAllSecond";
2988  // ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(TMath::Abs(finalMother->GetPdgCode()));
2989  // fillthis="particle_pdgAllSecondInvMass";
2990  // ((TH2F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(TMath::Abs(finalMother->GetPdgCode()),invariantMassMother);
2991  // if(finalMother->GetPdgCode()==511)
2992  // {
2993  // // isDesiredCandidate = kTRUE;
2994  // for (Int_t iDaughter = 0; iDaughter < finalMother->GetNDaughters(); ++iDaughter) //will work up to ten daughters
2995  // {
2996  // AliAODMCParticle* daughterfinalMother = (AliAODMCParticle*)mcTrackArray->At(finalMother->GetDaughter(0)+iDaughter);
2997  // fillthis="particle_daughterPdgTwoStep511";
2998  // ((TH2F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(iDaughter,TMath::Abs(daughterfinalMother->GetPdgCode()));
2999  // }
3000  // for (Int_t iDaughter = 0; iDaughter < D0GrandMother->GetNDaughters(); ++iDaughter) //will work up to ten daughters
3001  // {
3002  // AliAODMCParticle* daughterD0GrandMother = (AliAODMCParticle*)mcTrackArray->At(D0GrandMother->GetDaughter(0)+iDaughter);
3003  // fillthis="particle_daughterPdgTwoStep511";
3004  // ((TH2F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(iDaughter+10,TMath::Abs(daughterD0GrandMother->GetPdgCode()));
3005  // }
3006  // for (Int_t iDaughter = 0; iDaughter < D0Mother->GetNDaughters(); ++iDaughter) //will work up to ten daughters
3007  // {
3008  // AliAODMCParticle* daughterD0Mother = (AliAODMCParticle*)mcTrackArray->At(D0Mother->GetDaughter(0)+iDaughter);
3009  // fillthis="particle_daughterPdgTwoStep511";
3010  // ((TH2F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(iDaughter+20,TMath::Abs(daughterD0Mother->GetPdgCode()));
3011  // }
3012  // }
3013  // if(finalMother->GetPdgCode()==521)
3014  // {
3015  // // isDesiredCandidate = kTRUE;
3016  // for (Int_t iDaughter = 0; iDaughter < finalMother->GetNDaughters(); ++iDaughter) //will work up to ten daughters
3017  // {
3018  // AliAODMCParticle* daughterfinalMother = (AliAODMCParticle*)mcTrackArray->At(finalMother->GetDaughter(0)+iDaughter);
3019  // fillthis="particle_daughterPdgTwoStep521";
3020  // ((TH2F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(iDaughter,TMath::Abs(daughterfinalMother->GetPdgCode()));
3021  // }
3022  // for (Int_t iDaughter = 0; iDaughter < D0GrandMother->GetNDaughters(); ++iDaughter) //will work up to ten daughters
3023  // {
3024  // AliAODMCParticle* daughterD0GrandMother = (AliAODMCParticle*)mcTrackArray->At(D0GrandMother->GetDaughter(0)+iDaughter);
3025  // fillthis="particle_daughterPdgTwoStep521";
3026  // ((TH2F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(iDaughter+10,TMath::Abs(daughterD0GrandMother->GetPdgCode()));
3027  // }
3028  // for (Int_t iDaughter = 0; iDaughter < D0Mother->GetNDaughters(); ++iDaughter) //will work up to ten daughters
3029  // {
3030  // AliAODMCParticle* daughterD0Mother = (AliAODMCParticle*)mcTrackArray->At(D0Mother->GetDaughter(0)+iDaughter);
3031  // fillthis="particle_daughterPdgTwoStep521";
3032  // ((TH2F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(iDaughter+20,TMath::Abs(daughterD0Mother->GetPdgCode()));
3033  // }
3034  // }
3035  // }
3036  // }
3037  // }
3038 
3039  // }
3040  // }
3041 
3042 
3043  if (!bWrongSign)
3044  {
3045  histType = 2;
3046  FillBPlusHistograms(&trackBPlus, primaryVertex, bz, motherType, histType);
3047  if (fUseMCInfo && isDesiredCandidate)
3048  {
3049  //fill mc histograms
3050  histType = 3;
3051  FillBPlusHistograms(&trackBPlus, primaryVertex, bz, motherType, histType);
3052  }
3053  }
3054 
3055 
3056  // pdgMassMother=TDatabasePDG::Instance()->GetParticle(521)->Mass();
3057  // Double_t massWindow =fHistMassWindow; //GeV/c^2
3058  if (TMath::Abs(invariantMassMother - pdgMassMother) < massWindow)
3059  {
3060  if (!bWrongSign)
3061  {
3062  FillFinalTrackHistograms(&trackBPlus, primaryVertex, bz, isDesiredCandidate, mcTrackArray);
3063  if (!isDesiredCandidate)
3064  {
3065  motherType = 0; histType = 4; FillD0Histograms(trackD0, primaryVertex, bz, motherType, histType, pdgD0);
3066  motherType = 1; histType = 4; FillBPlusHistograms(&trackBPlus, primaryVertex, bz, motherType, histType);
3067  }
3068  if (isDesiredCandidate)
3069  {
3070  motherType = 0; histType = 5; FillD0Histograms(trackD0, primaryVertex, bz, motherType, histType, pdgD0);
3071  motherType = 1; histType = 5; FillBPlusHistograms(&trackBPlus, primaryVertex, bz, motherType, histType);
3072  }
3073  }
3074  }
3075 
3076 
3077  // Here we fill the histograms per pt bin and apply the same sign method
3078  TString ptBinMother = "";
3079  Int_t ptBin = fCuts->PtBin(trackBPlus.Pt());
3080  ptBinMother += "_ptbin_"; ptBinMother += fPtBinLimits[ptBin]; ptBinMother += "_to_"; ptBinMother += fPtBinLimits[ptBin+1];
3081  histType = 6 + 2 * ptBin;
3082 
3083  Int_t d0PtBin = fCuts->PtBinD0forD0ptbin(trackD0->Pt());
3084  Int_t histTypeD0 = 2 * d0PtBin;
3085 
3086 
3087  if (TMath::Abs(invariantMassMother - pdgMassMother) < massWindow)
3088  {
3089  if (!bWrongSign && histType > 5)
3090  {
3091  if (!isDesiredCandidate)
3092  {
3093  motherType = 0; FillD0Histograms(trackD0, primaryVertex, bz, motherType, histType, pdgD0);
3094  motherType = 1; FillBPlusHistograms(&trackBPlus, primaryVertex, bz, motherType, histType);
3095  motherType = 2; FillD0Histograms(trackD0, primaryVertex, bz, motherType, d0PtBin, pdgD0);
3096  }
3097 
3098  if (isDesiredCandidate)
3099  {
3100  histType += 1;
3101  motherType = 0; FillD0Histograms(trackD0, primaryVertex, bz, motherType, histType, pdgD0);
3102  motherType = 1; FillBPlusHistograms(&trackBPlus, primaryVertex, bz, motherType, histType);
3103  motherType = 2; FillD0Histograms(trackD0, primaryVertex, bz, motherType, d0PtBin + 1, pdgD0);
3104  }
3105  }
3106  }
3107 
3108  Double_t invmassDelta = DeltaInvMassBPlusKpipi(&trackBPlus);
3109  if(bWrongSign && iRot == 0)
3110  {
3111  fillthis="invariantMassBPlus";
3112  fillthis += "_SameSign";
3113  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3114  fillthis="invariantMassBPlus";
3115  fillthis += ptBinMother + "_SameSign";
3116  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3117  fillthis="invariantMassBPlus";
3118  fillthis += "_SignSum";
3119  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother,-1);
3120  fillthis="invariantMassBPlus";
3121  fillthis += ptBinMother + "_SignSum";
3122  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother,-1);
3123  }
3124  if(!bWrongSign)
3125  {
3126  if(iRot == 0)
3127  {
3128  fillthis="invariantMassBPlus";
3129  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3130  fillthis="invariantMassBPlus";
3131  fillthis += ptBinMother;
3132  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3133  fillthis="invariantMassBPlus";
3134  fillthis += "_SignSum";
3135  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother,1);
3136  fillthis="invariantMassBPlus";
3137  fillthis += ptBinMother + "_SignSum";
3138  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother,1);
3139 
3140  // fillthis = "invariantMassBPlusSignal_BA";
3141  // if(isDesiredCandidate) ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3142 
3143  // fillthis = "invariantMassBPlusCorrelated_BA";
3144  // if(bIsCorrelatedBackground) ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3145 
3146  // fillthis = "invariantMassBPlusBackground_BA";
3147  // if(!isDesiredCandidate && !bIsCorrelatedBackground) ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3148 
3149  // if(bIsCorrelatedBackground511)
3150  // {
3151  // fillthis="invariantMassBPlus_correlated511";
3152  // ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3153  // fillthis="invariantMassBPlus";
3154  // fillthis += ptBinMother + "_correlated511";
3155  // ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3156  // }
3157 
3158  if(!isDesiredCandidate && !bIsInjected)
3159  {
3160  TString signName = "_HIJING_Background";
3161  fillthis="invariantMassBPlus";
3162  fillthis += signName;
3163  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3164  fillthis="invariantMassBPlus";
3165  fillthis += ptBinMother + signName;
3166  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3167  }
3168  if(isDesiredCandidate && !bIsInjected)
3169  {
3170  TString signName = "_HIJING_Signal";
3171  fillthis="invariantMassBPlus";
3172  fillthis += signName;
3173  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3174  fillthis="invariantMassBPlus";
3175  fillthis += ptBinMother + signName;
3176  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3177  }
3178  }
3179  else
3180  {
3181  TString signName = "_Background_rotation";
3182  fillthis="invariantMassBPlus";
3183  fillthis += signName;
3184  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3185  fillthis="invariantMassBPlus";
3186  fillthis += ptBinMother + signName;
3187  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3188  if(!isDesiredCandidate && !bIsInjected)
3189  {
3190  signName = "_HIJING_Background_rotation";
3191  fillthis="invariantMassBPlus";
3192  fillthis += signName;
3193  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3194  fillthis="invariantMassBPlus";
3195  fillthis += ptBinMother + signName;
3196  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3197  }
3198  }
3199  }
3200 
3201  if(trackBPlus.Pt() > 6.0)
3202  {
3203  TString broadptBinMother = "_ptbin_6_to_inf";
3204  if(bWrongSign && iRot == 0)
3205  {
3206  fillthis="invariantMassBPlus";
3207  fillthis += broadptBinMother + "_SameSign";
3208  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3209  fillthis="invariantMassBPlus";
3210  fillthis += broadptBinMother + "_SignSum";
3211  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother,-1);
3212  }
3213  if(!bWrongSign)
3214  {
3215  if(iRot == 0)
3216  {
3217  fillthis="invariantMassBPlus";
3218  fillthis += broadptBinMother;
3219  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3220  fillthis="invariantMassBPlus";
3221  fillthis += broadptBinMother + "_SignSum";
3222  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother,1);
3223  if(!isDesiredCandidate && !bIsInjected)
3224  {
3225  TString signName = "_HIJING_Background";
3226  fillthis="invariantMassBPlus";
3227  fillthis += broadptBinMother + signName;
3228  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3229  }
3230  if(isDesiredCandidate && !bIsInjected)
3231  {
3232  TString signName = "_HIJING_Signal";
3233  fillthis="invariantMassBPlus";
3234  fillthis += broadptBinMother + signName;
3235  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3236  }
3237  }
3238  else
3239  {
3240  TString signName = "_Background_rotation";
3241  fillthis="invariantMassBPlus";
3242  fillthis += broadptBinMother + signName;
3243  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3244  if(!isDesiredCandidate && !bIsInjected)
3245  {
3246  signName = "_HIJING_Background_rotation";
3247  fillthis="invariantMassBPlus";
3248  fillthis += broadptBinMother + signName;
3249  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3250  }
3251  }
3252  }
3253  }
3254 
3255  if(trackBPlus.Pt() > 3.0)
3256  {
3257  TString broadptBinMother = "_ptbin_3_to_inf";
3258  if(bWrongSign && iRot == 0)
3259  {
3260  fillthis="invariantMassBPlus";
3261  fillthis += broadptBinMother + "_SameSign";
3262  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3263  fillthis="invariantMassBPlus";
3264  fillthis += broadptBinMother + "_SignSum";
3265  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother,-1);
3266  }
3267  if(!bWrongSign)
3268  {
3269  if(iRot == 0)
3270  {
3271  fillthis="invariantMassBPlus";
3272  fillthis += broadptBinMother;
3273  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3274  fillthis="invariantMassBPlus";
3275  fillthis += broadptBinMother + "_SignSum";
3276  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother,1);
3277  if(!isDesiredCandidate && !bIsInjected)
3278  {
3279  TString signName = "_HIJING_Background";
3280  fillthis="invariantMassBPlus";
3281  fillthis += broadptBinMother + signName;
3282  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3283  }
3284  if(isDesiredCandidate && !bIsInjected)
3285  {
3286  TString signName = "_HIJING_Signal";
3287  fillthis="invariantMassBPlus";
3288  fillthis += broadptBinMother + signName;
3289  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3290  }
3291  }
3292  else
3293  {
3294  TString signName = "_Background_rotation";
3295  fillthis="invariantMassBPlus";
3296  fillthis += broadptBinMother + signName;
3297  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3298  if(!isDesiredCandidate && !bIsInjected)
3299  {
3300  signName = "_HIJING_Background_rotation";
3301  fillthis="invariantMassBPlus";
3302  fillthis += broadptBinMother + signName;
3303  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invariantMassMother);
3304  }
3305  }
3306  }
3307  }
3308 
3309 
3310  //deltamass
3311  if(bWrongSign && iRot == 0)
3312  {
3313  fillthis="deltainvariantMassBPlus";
3314  fillthis += "_SameSign";
3315  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3316  fillthis="deltainvariantMassBPlus";
3317  fillthis += ptBinMother + "_SameSign";
3318  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3319  fillthis="deltainvariantMassBPlus";
3320  fillthis += "_SignSum";
3321  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta,-1);
3322  fillthis="deltainvariantMassBPlus";
3323  fillthis += ptBinMother + "_SignSum";
3324  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta,-1);
3325  }
3326  if(!bWrongSign)
3327  {
3328  if(iRot == 0)
3329  {
3330  fillthis="deltainvariantMassBPlus";
3331  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3332  fillthis="deltainvariantMassBPlus";
3333  fillthis += ptBinMother;
3334  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3335  fillthis="deltainvariantMassBPlus";
3336  fillthis += "_SignSum";
3337  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta,1);
3338  fillthis="deltainvariantMassBPlus";
3339  fillthis += ptBinMother + "_SignSum";
3340  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta,1);
3341 
3342  // if(bIsCorrelatedBackground511)
3343  // {
3344  // fillthis="deltainvariantMassBPlus_correlated511";
3345  // ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3346  // fillthis="deltainvariantMassBPlus";
3347  // fillthis += ptBinMother + "_correlated511" ;
3348  // ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3349  // }
3350 
3351  if(!isDesiredCandidate && !bIsInjected)
3352  {
3353  TString signName = "_HIJING_Background";
3354  fillthis="deltainvariantMassBPlus";
3355  fillthis += signName;
3356  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3357  fillthis="deltainvariantMassBPlus";
3358  fillthis += ptBinMother + signName;
3359  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3360  }
3361  if(isDesiredCandidate && !bIsInjected)
3362  {
3363  TString signName = "_HIJING_Signal";
3364  fillthis="deltainvariantMassBPlus";
3365  fillthis += signName;
3366  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3367  fillthis="deltainvariantMassBPlus";
3368  fillthis += ptBinMother + signName;
3369  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3370  }
3371  }
3372  else
3373  {
3374  TString signName = "_Background_rotation";
3375  fillthis="deltainvariantMassBPlus";
3376  fillthis += signName;
3377  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3378  fillthis="deltainvariantMassBPlus";
3379  fillthis += ptBinMother + signName;
3380  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3381  if(!isDesiredCandidate && !bIsInjected)
3382  {
3383  signName = "_HIJING_Background_rotation";
3384  fillthis="deltainvariantMassBPlus";
3385  fillthis += signName;
3386  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3387  fillthis="deltainvariantMassBPlus";
3388  fillthis += ptBinMother + signName;
3389  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3390  }
3391  }
3392  }
3393 
3394  if(trackBPlus.Pt() > 6.0)
3395  {
3396  TString broadptBinMother = "_ptbin_6_to_inf";
3397  if(bWrongSign && iRot == 0)
3398  {
3399  fillthis="deltainvariantMassBPlus";
3400  fillthis += broadptBinMother + "_SameSign";
3401  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3402  fillthis="deltainvariantMassBPlus";
3403  fillthis += broadptBinMother + "_SignSum";
3404  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta,-1);
3405  }
3406  if(!bWrongSign)
3407  {
3408  if(iRot == 0)
3409  {
3410  fillthis="deltainvariantMassBPlus";
3411  fillthis += broadptBinMother;
3412  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3413  fillthis="deltainvariantMassBPlus";
3414  fillthis += broadptBinMother + "_SignSum";
3415  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta,1);
3416  if(!isDesiredCandidate && !bIsInjected)
3417  {
3418  TString signName = "_HIJING_Background";
3419  fillthis="deltainvariantMassBPlus";
3420  fillthis += broadptBinMother + signName;
3421  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3422  }
3423  if(isDesiredCandidate && !bIsInjected)
3424  {
3425  TString signName = "_HIJING_Signal";
3426  fillthis="deltainvariantMassBPlus";
3427  fillthis += broadptBinMother + signName;
3428  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3429  }
3430  }
3431  else
3432  {
3433  TString signName = "_Background_rotation";
3434  fillthis="deltainvariantMassBPlus";
3435  fillthis += broadptBinMother + signName;
3436  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3437  if(!isDesiredCandidate && !bIsInjected)
3438  {
3439  signName = "_HIJING_Background_rotation";
3440  fillthis="deltainvariantMassBPlus";
3441  fillthis += broadptBinMother + signName;
3442  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3443  }
3444  }
3445  }
3446  }
3447 
3448  if(trackBPlus.Pt() > 3.0)
3449  {
3450  TString broadptBinMother = "_ptbin_3_to_inf";
3451  if(bWrongSign && iRot == 0)
3452  {
3453  fillthis="deltainvariantMassBPlus";
3454  fillthis += broadptBinMother + "_SameSign";
3455  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3456  fillthis="deltainvariantMassBPlus";
3457  fillthis += broadptBinMother + "_SignSum";
3458  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta,-1);
3459  }
3460  if(!bWrongSign)
3461  {
3462  if(iRot == 0)
3463  {
3464  fillthis="deltainvariantMassBPlus";
3465  fillthis += broadptBinMother;
3466  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3467  fillthis="deltainvariantMassBPlus";
3468  fillthis += broadptBinMother + "_SignSum";
3469  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta,1);
3470  if(!isDesiredCandidate && !bIsInjected)
3471  {
3472  TString signName = "_HIJING_Background";
3473  fillthis="deltainvariantMassBPlus";
3474  fillthis += broadptBinMother + signName;
3475  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3476  }
3477  if(isDesiredCandidate && !bIsInjected)
3478  {
3479  TString signName = "_HIJING_Signal";
3480  fillthis="deltainvariantMassBPlus";
3481  fillthis += broadptBinMother + signName;
3482  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3483  }
3484  }
3485  else
3486  {
3487  TString signName = "_Background_rotation";
3488  fillthis="deltainvariantMassBPlus";
3489  fillthis += broadptBinMother + signName;
3490  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3491  if(!isDesiredCandidate && !bIsInjected)
3492  {
3493  signName = "_HIJING_Background_rotation";
3494  fillthis="deltainvariantMassBPlus";
3495  fillthis += broadptBinMother + signName;
3496  ((TH1F*)(fOutputBPlusMC->FindObject(fillthis)))->Fill(invmassDelta);
3497  }
3498  }
3499  }
3500  }
3501 
3502  //we save the mother to a TClonesArray
3503  //if(!bWrongSign) new ((*fBPlusTracks)[iClonesArray++]) AliAODRecoDecayHF2Prong(*trackBPlus);
3504  delete vertexMother; vertexMother = NULL;
3505  delete trackBPlusPionRotated; trackBPlusPionRotated = nullptr;
3506  }
3507  }
3508 
3509  //Restore Object count
3510  //To save space in the table keeping track of all referenced objects,
3511  //we reset the object count to what it was at the beginning of the loop.
3512  TProcessID::SetObjectCount(ObjectNumber);
3513  }
3514  return;
3515 }
3516 //-------------------------------------------------------------------------------------
3517 void AliAnalysisTaskSEBPlustoD0Pi::FillFinalTrackHistograms(AliAODRecoDecayHF2Prong * selectedBPlus, AliAODVertex *primaryVertex, Double_t bz, Bool_t isDesiredCandidate,TClonesArray * mcTrackArray){
3518 
3519  //In this function we fill histograms with the properties of all the daughters of our selected signal candidate
3520 
3521  AliAODTrack* selectedBPlusPion = (AliAODTrack*)selectedBPlus->GetDaughter(0);
3522  AliAODRecoDecayHF2Prong* selectedD0 = (AliAODRecoDecayHF2Prong*)selectedBPlus->GetDaughter(1);
3523 
3524  AliAODTrack* selectedD0Pion;
3525  AliAODTrack* selectedD0Kaon;
3526 
3527  if(selectedBPlus->Charge() == 1) selectedD0Pion = (AliAODTrack*)selectedD0->GetDaughter(0);
3528  if(selectedBPlus->Charge() == -1) selectedD0Pion = (AliAODTrack*)selectedD0->GetDaughter(1);
3529 
3530  if(selectedBPlus->Charge() == 1) selectedD0Kaon = (AliAODTrack*)selectedD0->GetDaughter(1);
3531  if(selectedBPlus->Charge() == -1) selectedD0Kaon = (AliAODTrack*)selectedD0->GetDaughter(0);
3532 
3533  Double_t d0BPluspion = TMath::Abs(selectedBPlus->Getd0Prong(0));
3534  Double_t d0D0pion = 0;
3535  Double_t d0D0kaon = 0;
3536 
3537  if(selectedBPlus->Charge() == 1) d0D0pion = selectedD0->Getd0Prong(0);
3538  if(selectedBPlus->Charge() == -1) d0D0pion = selectedD0->Getd0Prong(1);
3539 
3540  if(selectedBPlus->Charge() == 1) d0D0kaon = selectedD0->Getd0Prong(1);
3541  if(selectedBPlus->Charge() == -1) d0D0kaon = selectedD0->Getd0Prong(0);
3542 
3543  Double_t pt_track = 0;
3544  Double_t momentum_track = 0;
3545  Int_t numberOfITS = 0;
3546  Int_t numberOfTPC = 0;
3547  Int_t daughterType, histType;
3548  Int_t totalNumberOfITS = 0;
3549  Int_t totalNumberOfTPC = 0;
3550  Double_t nSigmaTPC = 0;
3551  Double_t nSigmaTOF = 0;
3552  Double_t nSigmaTPCtotal = 0;
3553  Double_t nSigmaTOFtotal = 0;
3554  Int_t pionPIDnumber = 2;
3555  Int_t kaonPIDnumber = 3;
3556  Int_t TPCok = 0;
3557  Int_t TOFok = 0;
3558 
3559  AliAODPidHF* trackPIDHF = (AliAODPidHF*)fCuts->GetPidHF();
3560 
3561  //fill the D0 pion info
3562  pt_track = selectedD0Pion->Pt();
3563  momentum_track = selectedD0Pion->P();
3564  numberOfITS = selectedD0Pion->GetITSNcls();
3565  numberOfTPC = selectedD0Pion->GetTPCNcls();
3566  totalNumberOfITS += numberOfITS;
3567  totalNumberOfTPC += numberOfTPC;
3568  TPCok = trackPIDHF->GetnSigmaTPC(selectedD0Pion, pionPIDnumber, nSigmaTPC);
3569  TOFok = trackPIDHF->GetnSigmaTOF(selectedD0Pion, pionPIDnumber, nSigmaTOF);
3570  if(TPCok != -1) nSigmaTPCtotal += nSigmaTPC*nSigmaTPC;
3571  if(TOFok != -1) nSigmaTOFtotal += nSigmaTOF*nSigmaTOF;
3572 
3573  Double_t ptBPlus = selectedBPlus->Pt();
3574 
3575  daughterType = 0;
3576  histType = 4;
3577  if(!isDesiredCandidate)
3578  {
3579  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
3580  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
3581  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
3582  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
3583 
3584  for (Int_t j = 0; j < 10; ++j)
3585  {
3586  if(selectedD0Pion->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
3587 
3588  }
3589 
3590  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
3591  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
3592  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
3593  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0D0pion);
3594  ((TH1F*)fDaughterHistogramArray[daughterType][histType][9])->Fill(selectedD0Pion->Eta());
3595  ((TH2F*)fDaughterHistogramArray2D[daughterType][4])->Fill(ptBPlus,pt_track);
3596  }
3597 
3598  if(isDesiredCandidate)
3599  {
3600  histType = 5;
3601  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
3602  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
3603  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
3604  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
3605 
3606  for (Int_t j = 0; j < 10; ++j)
3607  {
3608  if(selectedD0Pion->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
3609 
3610  }
3611 
3612  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
3613  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
3614  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
3615  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0D0pion);
3616  ((TH1F*)fDaughterHistogramArray[daughterType][histType][9])->Fill(selectedD0Pion->Eta());
3617  ((TH2F*)fDaughterHistogramArray2D[daughterType][5])->Fill(ptBPlus,pt_track);
3618  }
3619 
3620  //we save the pdgcode of the used particle and its mother to check PID efficiency
3621  if(fUseMCInfo)
3622  {
3623  Float_t pdgCodeParticle = -1;
3624  Float_t pdgCodeParticleMother = -1;
3625  Int_t mcLabelParticle = -1;
3626  Int_t mcLabelParticleMother = -1;
3627  mcLabelParticle = selectedD0Pion->GetLabel();
3628 
3629  if(mcLabelParticle >= 0){
3630 
3631  AliAODMCParticle *mcTrackParticle = (AliAODMCParticle*)mcTrackArray->At(mcLabelParticle);
3632  pdgCodeParticle = TMath::Abs(mcTrackParticle->GetPdgCode());
3633  ((TH1F*)fDaughterHistogramArrayExtra[0][2])->Fill(pdgCodeParticle);
3634  mcLabelParticleMother = mcTrackParticle->GetMother();
3635 
3636  if(mcLabelParticleMother >= 0){
3637  AliAODMCParticle *mcTrackParticleMother = (AliAODMCParticle*)mcTrackArray->At(mcLabelParticleMother);
3638  pdgCodeParticleMother = TMath::Abs(mcTrackParticleMother->GetPdgCode());
3639  ((TH1F*)fDaughterHistogramArrayExtra[0][3])->Fill(pdgCodeParticleMother);
3640  }
3641  }
3642  }
3643 
3644  //fill the D0 kaon info
3645  pt_track = selectedD0Kaon->Pt();
3646  momentum_track = selectedD0Kaon->P();
3647  numberOfITS = selectedD0Kaon->GetITSNcls();
3648  numberOfTPC = selectedD0Kaon->GetTPCNcls();
3649  totalNumberOfITS += numberOfITS;
3650  totalNumberOfTPC += numberOfTPC;
3651  TPCok = trackPIDHF->GetnSigmaTPC(selectedD0Kaon, kaonPIDnumber, nSigmaTPC);
3652  TOFok = trackPIDHF->GetnSigmaTOF(selectedD0Kaon, kaonPIDnumber, nSigmaTOF);
3653  if(TPCok != -1) nSigmaTPCtotal += nSigmaTPC*nSigmaTPC;
3654  if(TOFok != -1) nSigmaTOFtotal += nSigmaTOF*nSigmaTOF;
3655 
3656  daughterType = 1;
3657  histType = 4;
3658  if(!isDesiredCandidate)
3659  {
3660  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
3661  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
3662  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
3663  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
3664 
3665  for (Int_t j = 0; j < 10; ++j)
3666  {
3667  if(selectedD0Kaon->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
3668 
3669  }
3670 
3671  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
3672  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
3673  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
3674  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0D0kaon);
3675  ((TH1F*)fDaughterHistogramArray[daughterType][histType][9])->Fill(selectedD0Kaon->Eta());
3676  ((TH2F*)fDaughterHistogramArray2D[daughterType][4])->Fill(ptBPlus,pt_track);
3677  }
3678 
3679  if(isDesiredCandidate)
3680  {
3681  histType = 5;
3682  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
3683  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
3684  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
3685  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
3686 
3687  for (Int_t j = 0; j < 10; ++j)
3688  {
3689  if(selectedD0Kaon->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
3690 
3691  }
3692 
3693  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
3694  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
3695  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
3696  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0D0kaon);
3697  ((TH1F*)fDaughterHistogramArray[daughterType][histType][9])->Fill(selectedD0Kaon->Eta());
3698  ((TH2F*)fDaughterHistogramArray2D[daughterType][5])->Fill(ptBPlus,pt_track);
3699  }
3700 
3701  //we save the pdgcode of the used particle and its mother to check PID efficiency
3702  if(fUseMCInfo)
3703  {
3704  Float_t pdgCodeParticle = -1;
3705  Float_t pdgCodeParticleMother = -1;
3706  Int_t mcLabelParticle = -1;
3707  Int_t mcLabelParticleMother = -1;
3708  mcLabelParticle = selectedD0Kaon->GetLabel();
3709 
3710  if(mcLabelParticle >= 0){
3711 
3712  AliAODMCParticle *mcTrackParticle = (AliAODMCParticle*)mcTrackArray->At(mcLabelParticle);
3713  pdgCodeParticle = TMath::Abs(mcTrackParticle->GetPdgCode());
3714  ((TH1F*)fDaughterHistogramArrayExtra[1][2])->Fill(pdgCodeParticle);
3715  mcLabelParticleMother = mcTrackParticle->GetMother();
3716 
3717  if(mcLabelParticleMother >= 0){
3718  AliAODMCParticle *mcTrackParticleMother = (AliAODMCParticle*)mcTrackArray->At(mcLabelParticleMother);
3719  pdgCodeParticleMother = TMath::Abs(mcTrackParticleMother->GetPdgCode());
3720  ((TH1F*)fDaughterHistogramArrayExtra[1][3])->Fill(pdgCodeParticleMother);
3721  }
3722  }
3723  }
3724 
3725  //fill the BPlus pion info
3726  pt_track = selectedBPlusPion->Pt();
3727  momentum_track = selectedBPlusPion->P();
3728  numberOfITS = selectedBPlusPion->GetITSNcls();
3729  numberOfTPC = selectedBPlusPion->GetTPCNcls();
3730  totalNumberOfITS += numberOfITS;
3731  totalNumberOfTPC += numberOfTPC;
3732  TPCok = trackPIDHF->GetnSigmaTPC(selectedBPlusPion, pionPIDnumber, nSigmaTPC);
3733  TOFok = trackPIDHF->GetnSigmaTOF(selectedBPlusPion, pionPIDnumber, nSigmaTOF);
3734  if(TPCok != -1) nSigmaTPCtotal += nSigmaTPC*nSigmaTPC;
3735  if(TOFok != -1) nSigmaTOFtotal += nSigmaTOF*nSigmaTOF;
3736 
3737  daughterType = 2;
3738  histType = 4;
3739  if(!isDesiredCandidate)
3740  {
3741  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
3742  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
3743  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
3744  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
3745 
3746  for (Int_t j = 0; j < 10; ++j)
3747  {
3748  if(selectedBPlusPion->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
3749 
3750  }
3751 
3752  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
3753  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
3754  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
3755  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0BPluspion);
3756  ((TH1F*)fDaughterHistogramArray[daughterType][histType][9])->Fill(selectedBPlusPion->Eta());
3757  ((TH2F*)fDaughterHistogramArray2D[daughterType][4])->Fill(ptBPlus,pt_track);
3758  }
3759 
3760  if(isDesiredCandidate)
3761  {
3762  histType = 5;
3763  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
3764  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
3765  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
3766  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
3767 
3768  for (Int_t j = 0; j < 10; ++j)
3769  {
3770  if(selectedBPlusPion->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
3771 
3772  }
3773 
3774  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
3775  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
3776  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
3777  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0BPluspion);
3778  ((TH1F*)fDaughterHistogramArray[daughterType][histType][9])->Fill(selectedBPlusPion->Eta());
3779  ((TH2F*)fDaughterHistogramArray2D[daughterType][5])->Fill(ptBPlus,pt_track);
3780  }
3781 
3782  //we save the pdgcode of the used particle and its mother to check PID efficiency
3783  if(fUseMCInfo)
3784  {
3785  Float_t pdgCodeParticle = -1;
3786  Float_t pdgCodeParticleMother = -1;
3787  Int_t mcLabelParticle = -1;
3788  Int_t mcLabelParticleMother = -1;
3789  mcLabelParticle = selectedBPlusPion->GetLabel();
3790 
3791  if(mcLabelParticle >= 0){
3792 
3793  AliAODMCParticle *mcTrackParticle = (AliAODMCParticle*)mcTrackArray->At(mcLabelParticle);
3794  pdgCodeParticle = TMath::Abs(mcTrackParticle->GetPdgCode());
3795  ((TH1F*)fDaughterHistogramArrayExtra[daughterType][2])->Fill(pdgCodeParticle);
3796  mcLabelParticleMother = mcTrackParticle->GetMother();
3797 
3798  if(mcLabelParticleMother >= 0){
3799  AliAODMCParticle *mcTrackParticleMother = (AliAODMCParticle*)mcTrackArray->At(mcLabelParticleMother);
3800  pdgCodeParticleMother = TMath::Abs(mcTrackParticleMother->GetPdgCode());
3801  ((TH1F*)fDaughterHistogramArrayExtra[daughterType][3])->Fill(pdgCodeParticleMother);
3802  }
3803  }
3804  }
3805 
3806  if(!isDesiredCandidate)
3807  {
3808  ((TH1F*)(fOutputBPlusMC->FindObject("totalITSBackground")))->Fill(totalNumberOfITS);
3809  ((TH1F*)(fOutputBPlusMC->FindObject("totalTPCBackground")))->Fill(totalNumberOfTPC);
3810  ((TH1F*)(fOutputBPlusMC->FindObject("totalSigmaPIDBackground")))->Fill(sqrt(nSigmaTPCtotal + nSigmaTOFtotal));
3811  }
3812  if(isDesiredCandidate)
3813  {
3814  ((TH1F*)(fOutputBPlusMC->FindObject("totalITSSignal")))->Fill(totalNumberOfITS);
3815  ((TH1F*)(fOutputBPlusMC->FindObject("totalTPCSignal")))->Fill(totalNumberOfTPC);
3816  ((TH1F*)(fOutputBPlusMC->FindObject("totalSigmaPIDSignal")))->Fill(sqrt(nSigmaTPCtotal + nSigmaTOFtotal));
3817  }
3818 
3819  // we look at the invariant mass combinations of all daughter tracks
3820  AliExternalTrackParam trackBPlusPion;
3821  trackBPlusPion.CopyFromVTrack(selectedBPlusPion);
3822  AliExternalTrackParam trackD0Pion;
3823  trackD0Pion.CopyFromVTrack(selectedD0Pion);
3824  AliExternalTrackParam trackD0Kaon;
3825  trackD0Kaon.CopyFromVTrack(selectedD0Kaon);
3826 
3827  UInt_t prongs2[2] = {0};
3828  UInt_t prongs3[3] = {0};
3829 
3830  prongs2[0] = 211; prongs2[1] = 211;
3831  TwoTrackCombinationInfo(&trackD0Pion, &trackBPlusPion, primaryVertex, bz, isDesiredCandidate, "invmassD0PionBPlusPion", prongs2);
3832  prongs2[0] = 321; prongs2[1] = 211;
3833  TwoTrackCombinationInfo(&trackD0Kaon, &trackBPlusPion, primaryVertex, bz, isDesiredCandidate, "invmassD0KaonBPlusPion", prongs2);
3834  prongs3[0] = 211; prongs3[1] = 321; prongs3[2] = 211;
3835  ThreeTrackCombinationInfo(&trackD0Pion, &trackD0Kaon, &trackBPlusPion, primaryVertex, bz, isDesiredCandidate, "invmassD0PionD0KaonBPlusPion", prongs3);
3836 
3837  return;
3838 }
3839 //-------------------------------------------------------------------------------------
3841 {
3845 
3846  AliAODRecoDecayHF2Prong * D0 = (AliAODRecoDecayHF2Prong*)BPlus->GetDaughter(1);
3847 
3848  Double_t e[3] = {0};
3849  if (BPlus->Charge() == -1)
3850  {
3851  e[0] = D0->EProng(0, 211);
3852  e[1] = D0->EProng(1, 321);
3853  } else if (BPlus->Charge() == 1) {
3854  e[0] = D0->EProng(0, 321);
3855  e[1] = D0->EProng(1, 211);
3856  }
3857  e[2] = BPlus->EProng(0, 211);
3858 
3859  Double_t esum = e[0] + e[1] + e[2];
3860  Double_t invMassBPlus = TMath::Sqrt(esum * esum - BPlus->P2());
3861 
3862  Double_t invMassD0 = -1;
3863 
3864  if (BPlus->Charge() == -1) {invMassD0 = D0->InvMassD0();}
3865  else {invMassD0 = D0->InvMassD0bar();}
3866  if (invMassD0 == -1) {std::cout << "wrong invmass delta D0 BPlus" << std::endl;}
3867 
3868  return invMassBPlus - invMassD0;
3869 }
3870 //-------------------------------------------------------------------------------------
3871 void AliAnalysisTaskSEBPlustoD0Pi::FillD0Histograms(AliAODRecoDecayHF2Prong * selectedMother, AliAODVertex *primaryVertex, Double_t bz, Int_t motherType, Int_t histType, Int_t pdgCodeMother) {
3872 
3873  if(histType<0) return;
3874 
3875  //In this function we fill the histograms of the reconstructed mothers
3876  Double_t ptMother = 0.0;
3877  Double_t momentumMother = 0.0;
3878  Double_t etaMother = 0.0;
3879  Double_t phiMother = 0.0;
3880  Double_t d0Mother = 0.0;
3881  Double_t d0firstTrack = 0.0;
3882  Double_t d0secondTrack = 0.0;
3883  Double_t pointingAngle = 0.0;
3884  Double_t impactProduct = 0.0;
3885  Double_t impactProductXY = 0.0;
3886  Double_t invariantMassMother = 0.0;
3887  Double_t invmassDelta = 0.0;
3888  Double_t dcaMother = 0.0;
3889  AliAODVertex * vertexMother = 0x0;
3890  AliAODVertex * vertexDaughter = 0x0;
3891  Double_t vertexDistance = 0.0;
3892  Double_t decayLengthDaughter = 0.0;
3893  Double_t decayTime = 0.0;
3894  Double_t angleMotherFirstDaughter = 0.0;
3895  Double_t angleMotherSecondDaughter = 0.0;
3896  Double_t topomaticFirstDaughter = 0.0;
3897  Double_t topomaticSecondDaughter = 0.0;
3898  Double_t ptFirstDaughter = 0.0;
3899  Double_t ptSecondDaughter = 0.0;
3900  UInt_t prongs[2], prongs2[2];
3901  Double_t angleBetweenBothDaughters = 0;
3902  Double_t cosThetaStar = 0;
3903  Double_t normDecayLength = 0;
3904  Double_t pdgMassMother=TDatabasePDG::Instance()->GetParticle(421)->Mass();
3905 
3906 
3907  prongs[0] = 211; prongs[1] = 321;
3908  prongs2[0] = 321; prongs2[1] = 211;
3909  AliAODTrack * firstDaughter = (AliAODTrack*)selectedMother->GetDaughter(0);
3910  AliAODTrack * secondDaughter = (AliAODTrack*)selectedMother->GetDaughter(1);
3911  vertexMother = selectedMother->GetSecondaryVtx();
3912  ptFirstDaughter = firstDaughter->Pt();
3913  ptSecondDaughter = secondDaughter->Pt();
3914 
3915  //Topomatic
3916  Double_t dd0pr1=0.;
3917  Double_t dd0pr2=0.;
3918  Double_t dd0max=0.;
3919  Double_t dd0min=0.;
3920  for(Int_t ipr=0; ipr<2; ipr++)
3921  {
3922  Double_t diffIP, errdiffIP;
3923  selectedMother->Getd0MeasMinusExpProng(ipr,bz,diffIP,errdiffIP);
3924  Double_t normdd0=0.;
3925  if(errdiffIP>0.) normdd0=diffIP/errdiffIP;
3926  if(ipr==0) dd0pr1=normdd0;
3927  if(ipr==1) dd0pr2=normdd0;
3928 
3929  }
3930  if(TMath::Abs(dd0pr1)>TMath::Abs(dd0pr2)) {dd0max=dd0pr1; dd0min=dd0pr2;}
3931  else {dd0max=dd0pr2; dd0min=dd0pr1;}
3932 
3933  AliExternalTrackParam motherTrack;
3934  motherTrack.CopyFromVTrack(selectedMother);
3935  Double_t d0z0[2],covd0z0[3],d0[2];
3936  motherTrack.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
3937  d0[0] = d0z0[0];
3938 
3939  ptMother = selectedMother->Pt();
3940  momentumMother = selectedMother->P();
3941  etaMother = selectedMother->Eta();
3942  phiMother = selectedMother->Phi();
3943 
3944  d0Mother = TMath::Abs(d0[0]);
3945  d0firstTrack = TMath::Abs(selectedMother->Getd0Prong(0));
3946  d0secondTrack = TMath::Abs(selectedMother->Getd0Prong(1));
3947  pointingAngle = selectedMother->CosPointingAngle();
3948  impactProduct = selectedMother->Prodd0d0();
3949  impactProductXY = TMath::Abs(selectedMother->ImpParXY());
3950  invariantMassMother = selectedMother->InvMass(2,prongs);
3951  if(pdgCodeMother == -421) invariantMassMother = selectedMother->InvMass(2,prongs2);
3952 
3953  dcaMother = selectedMother->GetDCA();
3954  vertexDistance = vertexMother->DistanceToVertex(primaryVertex);
3955  angleMotherFirstDaughter = (selectedMother->Px() * firstDaughter->Px() + selectedMother->Py() * firstDaughter->Py() + selectedMother->Pz() * firstDaughter->Pz()) /(selectedMother->P() * firstDaughter->P());
3956  angleMotherSecondDaughter = (selectedMother->Px() * secondDaughter->Px() + selectedMother->Py() * secondDaughter->Py() + selectedMother->Pz() * secondDaughter->Pz()) /(selectedMother->P() * secondDaughter->P());
3957  cosThetaStar = selectedMother->CosThetaStar(0,421,211,321);
3958  angleBetweenBothDaughters = (firstDaughter->Px() * secondDaughter->Px() + firstDaughter->Py() * secondDaughter->Py() + firstDaughter->Pz() * secondDaughter->Pz()) /(firstDaughter->P() * secondDaughter->P());
3959  normDecayLength = selectedMother->NormalizedDecayLength();
3960 
3961  Double_t pseudoProperDecayLength = ((vertexMother->GetX() - primaryVertex->GetX()) * selectedMother->Px() / TMath::Abs(selectedMother->Pt())) + ((vertexMother->GetY() - primaryVertex->GetY()) * selectedMother->Py() / TMath::Abs(selectedMother->Pt()));
3962  Double_t pseudoProperDecayTime = pseudoProperDecayLength * pdgMassMother/ptMother;
3963  decayTime = vertexDistance / (299792458 * TMath::Sqrt(1/((pdgMassMother*pdgMassMother/(momentumMother*momentumMother)) + 1)));
3964 
3965  Double_t phi = selectedMother->Phi();
3966  Double_t theta = selectedMother->Theta();
3967  Double_t covMatrix[21];
3968  selectedMother->GetCovarianceXYZPxPyPz(covMatrix);
3969 
3970  Double_t cp = TMath::Cos(phi);
3971  Double_t sp = TMath::Sin(phi);
3972  Double_t ct = TMath::Cos(theta);
3973  Double_t st = TMath::Sin(theta);
3974 
3975  Double_t errorMomentum = covMatrix[9]*cp*cp*ct*ct // GetCovPxPx
3976  +covMatrix[13]*2.*cp*sp*ct*ct // GetCovPxPy
3977  +covMatrix[18]*2.*cp*ct*st // GetCovPxPz
3978  +covMatrix[14]*sp*sp*ct*ct // GetCovPyPy
3979  +covMatrix[19]*2.*sp*ct*st // GetCovPyPz
3980  +covMatrix[20]*st*st; // GetCovPzPz
3981  Double_t normalizedDecayTime = selectedMother->NormalizedDecayLength() / (299792458 * TMath::Sqrt(1/((pdgMassMother*pdgMassMother*errorMomentum*errorMomentum/(momentumMother*momentumMother)) + 1)));
3982 
3983  Double_t eKaon = selectedMother->EProng(1,321);
3984  Double_t invMassKaon = TMath::Sqrt(eKaon*eKaon-secondDaughter->P()*secondDaughter->P());
3985  Double_t invMassD0 = selectedMother->InvMassD0();
3986  invmassDelta = invMassD0 - invMassKaon;
3987 
3988  Double_t vertexMotherX = vertexMother->GetX();
3989  Double_t vertexMotherY = vertexMother->GetY();
3990  Double_t vertexMotherZ = vertexMother->GetZ();
3991 
3992  Double_t cosPointingAngleXY = selectedMother->CosPointingAngleXY();
3993  Double_t distanceXYToVertex = vertexMother->DistanceXYToVertex(primaryVertex);
3994  Double_t normalizedDecayLengthXY = selectedMother->NormalizedDecayLengthXY();
3995 
3996  ((TH1F*)fMotherHistogramArray[motherType][histType][0])->Fill(ptMother);
3997  ((TH1F*)fMotherHistogramArray[motherType][histType][1])->Fill(ptFirstDaughter);
3998  ((TH1F*)fMotherHistogramArray[motherType][histType][2])->Fill(ptSecondDaughter);
3999  ((TH1F*)fMotherHistogramArray[motherType][histType][3])->Fill(etaMother);
4000  ((TH1F*)fMotherHistogramArray[motherType][histType][4])->Fill(phiMother);
4001  ((TH1F*)fMotherHistogramArray[motherType][histType][5])->Fill(d0Mother);
4002  ((TH1F*)fMotherHistogramArray[motherType][histType][6])->Fill(d0firstTrack);
4003  ((TH1F*)fMotherHistogramArray[motherType][histType][7])->Fill(d0secondTrack);
4004  ((TH1F*)fMotherHistogramArray[motherType][histType][8])->Fill(pointingAngle);
4005  ((TH1F*)fMotherHistogramArray[motherType][histType][9])->Fill(impactProduct);
4006  ((TH1F*)fMotherHistogramArray[motherType][histType][10])->Fill(impactProductXY);
4007  ((TH1F*)fMotherHistogramArray[motherType][histType][11])->Fill(invariantMassMother);
4008  ((TH1F*)fMotherHistogramArray[motherType][histType][12])->Fill(invmassDelta);
4009  ((TH1F*)fMotherHistogramArray[motherType][histType][13])->Fill(dcaMother);
4010  ((TH1F*)fMotherHistogramArray[motherType][histType][14])->Fill(vertexDistance);
4011  ((TH1F*)fMotherHistogramArray[motherType][histType][15])->Fill(normDecayLength);
4012  ((TH1F*)fMotherHistogramArray[motherType][histType][16])->Fill(pseudoProperDecayTime);
4013  ((TH1F*)fMotherHistogramArray[motherType][histType][17])->Fill(decayTime);
4014  ((TH1F*)fMotherHistogramArray[motherType][histType][18])->Fill(normalizedDecayTime);
4015  ((TH1F*)fMotherHistogramArray[motherType][histType][19])->Fill(angleMotherFirstDaughter);
4016  ((TH1F*)fMotherHistogramArray[motherType][histType][20])->Fill(angleMotherSecondDaughter);
4017  ((TH1F*)fMotherHistogramArray[motherType][histType][21])->Fill(angleBetweenBothDaughters);
4018  ((TH1F*)fMotherHistogramArray[motherType][histType][22])->Fill(cosThetaStar);
4019 
4020  ((TH1F*)fMotherHistogramArray[motherType][histType][23])->Fill(vertexMotherX);
4021  ((TH1F*)fMotherHistogramArray[motherType][histType][24])->Fill(vertexMotherY);
4022  ((TH1F*)fMotherHistogramArray[motherType][histType][25])->Fill(vertexMotherZ);
4023 
4024  ((TH1F*)fMotherHistogramArray[motherType][histType][36])->Fill(TMath::Abs(dd0pr1));
4025  ((TH1F*)fMotherHistogramArray[motherType][histType][37])->Fill(TMath::Abs(dd0pr2));
4026  ((TH1F*)fMotherHistogramArray[motherType][histType][38])->Fill(TMath::Abs(dd0max));
4027  ((TH1F*)fMotherHistogramArray[motherType][histType][39])->Fill(TMath::Abs(dd0min));
4028 
4029  ((TH1F*)fMotherHistogramArray[motherType][histType][40])->Fill(cosPointingAngleXY);
4030  ((TH1F*)fMotherHistogramArray[motherType][histType][41])->Fill(distanceXYToVertex);
4031  ((TH1F*)fMotherHistogramArray[motherType][histType][42])->Fill(normalizedDecayLengthXY);
4032  ((TH1F*)fMotherHistogramArray[motherType][histType][43])->Fill(vertexMother->GetChi2());
4033  ((TH1F*)fMotherHistogramArray[motherType][histType][44])->Fill(vertexMother->GetChi2perNDF());
4034 
4035 
4036  //we fill the 2D histograms
4037  Int_t nFirst = 0;
4038  Int_t nSecond = 1;
4039  Int_t nVariables = 10;
4040  Int_t nHistograms = nVariables * (nVariables - 1) / 2;
4041  for (Int_t k = 0; k < nHistograms; ++k)
4042  {
4043  Double_t firstVariable = 0.0;
4044  Double_t secondVariable = 0.0;
4045 
4046  if(nFirst==0) firstVariable = d0firstTrack;
4047  if(nFirst==1) firstVariable = d0secondTrack;
4048  if(nFirst==2) firstVariable = d0Mother;
4049  if(nFirst==3) firstVariable = pointingAngle;
4050  if(nFirst==4) firstVariable = impactProduct;
4051  if(nFirst==5) firstVariable = impactProductXY;
4052  if(nFirst==6) firstVariable = vertexDistance;
4053  if(nFirst==7) firstVariable = normDecayLength;
4054  if(nFirst==8) firstVariable = cosPointingAngleXY;
4055  if(nFirst==9) firstVariable = distanceXYToVertex;
4056  if(nFirst==10) firstVariable = normalizedDecayLengthXY;
4057 
4058  if(nSecond==0) secondVariable = d0firstTrack;
4059  if(nSecond==1) secondVariable = d0secondTrack;
4060  if(nSecond==2) secondVariable = d0Mother;
4061  if(nSecond==3) secondVariable = pointingAngle;
4062  if(nSecond==4) secondVariable = impactProduct;
4063  if(nSecond==5) secondVariable = impactProductXY;
4064  if(nSecond==6) secondVariable = vertexDistance;
4065  if(nSecond==7) secondVariable = normDecayLength;
4066  if(nSecond==8) secondVariable = cosPointingAngleXY;
4067  if(nSecond==9) secondVariable = distanceXYToVertex;
4068  if(nSecond==10) secondVariable = normalizedDecayLengthXY;
4069 
4070  ((TH2F*)fMotherHistogramArray2D[motherType][histType][k])->Fill(firstVariable,secondVariable);
4071 
4072  nSecond++;
4073  if(nSecond>nVariables)
4074  {
4075  nFirst++;
4076  nSecond = nFirst + 1;
4077  }
4078  }
4079 
4080  return;
4081 }
4082 //-------------------------------------------------------------------------------------
4083 void AliAnalysisTaskSEBPlustoD0Pi::FillBPlusHistograms(AliAODRecoDecayHF2Prong * selectedMother, AliAODVertex *primaryVertex, Double_t bz, Int_t motherType, Int_t histType) {
4084 
4085  //In this function we fill the histograms of the reconstructed mothers
4086  Double_t ptMother = 0.0;
4087  Double_t momentumMother = 0.0;
4088  Double_t etaMother = 0.0;
4089  Double_t phiMother = 0.0;
4090  Double_t d0Mother = 0.0;
4091  Double_t d0firstTrack = 0.0;
4092  Double_t d0secondTrack = 0.0;
4093  Double_t pointingAngle = 0.0;
4094  Double_t impactProduct = 0.0;
4095  Double_t impactProductXY = 0.0;
4096  Double_t invariantMassMother = 0.0;
4097  Double_t invmassDelta = 0.0;
4098  Double_t dcaMother = 0.0;
4099  AliAODVertex * vertexMother = 0x0;
4100  AliAODVertex * vertexDaughter = 0x0;
4101  Double_t vertexDistance = 0.0;
4102  Double_t normDecayLength = 0.0;
4103  Double_t decayTime = 0.0;
4104  Double_t angleMotherFirstDaughter = 0.0;
4105  Double_t angleMotherSecondDaughter = 0.0;
4106  Double_t topomaticFirstDaughter = 0.0;
4107  Double_t topomaticSecondDaughter = 0.0;
4108  Double_t ptFirstDaughter = 0.0;
4109  Double_t ptSecondDaughter = 0.0;
4110  UInt_t prongs[2];
4111  Double_t cosThetaStar = 0;
4112  Double_t angleBetweenBothDaughters = 0;
4113  Double_t d0MultiProduct = 0;
4114 
4115  Double_t pdgMassMother = 0;
4116 
4117  AliExternalTrackParam motherTrack;
4118  motherTrack.CopyFromVTrack(selectedMother);
4119  Double_t d0z0[2], covd0z0[3], d0[2];
4120  motherTrack.PropagateToDCA(primaryVertex, bz, 100., d0z0, covd0z0);
4121  d0[0] = d0z0[0];
4122 
4123  //BPlus
4124  prongs[1] = 421; prongs[0] = 211;
4125  invmassDelta = DeltaInvMassBPlusKpipi(selectedMother);
4126  AliAODTrack * firstDaughter = (AliAODTrack*)selectedMother->GetDaughter(0);
4127  AliAODRecoDecayHF2Prong * secondDaughter = (AliAODRecoDecayHF2Prong*)selectedMother->GetDaughter(1);
4128 
4129  ptFirstDaughter = firstDaughter->Pt();
4130  ptSecondDaughter = secondDaughter->Pt();
4131  vertexMother = selectedMother->GetSecondaryVtx();
4132  vertexDaughter = secondDaughter->GetSecondaryVtx();
4133  angleMotherFirstDaughter = (selectedMother->Px() * firstDaughter->Px() + selectedMother->Py() * firstDaughter->Py() + selectedMother->Pz() * firstDaughter->Pz()) / (selectedMother->P() * firstDaughter->P());
4134  angleMotherSecondDaughter = (selectedMother->Px() * secondDaughter->Px() + selectedMother->Py() * secondDaughter->Py() + selectedMother->Pz() * secondDaughter->Pz()) / (selectedMother->P() * secondDaughter->P());
4135  angleBetweenBothDaughters = (firstDaughter->Px() * secondDaughter->Px() + firstDaughter->Py() * secondDaughter->Py() + firstDaughter->Pz() * secondDaughter->Pz()) / (firstDaughter->P() * secondDaughter->P());
4136  cosThetaStar = selectedMother->CosThetaStar(0, 521, 421, 211);
4137  pdgMassMother = TDatabasePDG::Instance()->GetParticle(521)->Mass();
4138  d0MultiProduct = selectedMother->Getd0Prong(0) * secondDaughter->Getd0Prong(0) * secondDaughter->Getd0Prong(1);
4139 
4140  //Topomatic
4141  Double_t dd0pr1 = 0.;
4142  Double_t dd0pr2 = 0.;
4143  Double_t dd0max = 0.;
4144  Double_t dd0min = 0.;
4145  for (Int_t ipr = 0; ipr < 2; ipr++)
4146  {
4147  Double_t diffIP, errdiffIP;
4148  selectedMother->Getd0MeasMinusExpProng(ipr, bz, diffIP, errdiffIP);
4149  Double_t normdd0 = 0.;
4150  if (errdiffIP > 0.) normdd0 = diffIP / errdiffIP;
4151  if (ipr == 0) dd0pr1 = normdd0;
4152  if (ipr == 1) dd0pr2 = normdd0;
4153 
4154  // else if(TMath::Abs(normdd0)>TMath::Abs(dd0max)) dd0max=normdd0;
4155  }
4156  if (TMath::Abs(dd0pr1) > TMath::Abs(dd0pr2)) {dd0max = dd0pr1; dd0min = dd0pr2;}
4157  else {dd0max = dd0pr2; dd0min = dd0pr1;}
4158 
4159 
4160 
4161  ptMother = selectedMother->Pt();
4162  momentumMother = selectedMother->P();
4163  etaMother = selectedMother->Eta();
4164  phiMother = selectedMother->Phi();
4165  d0Mother = TMath::Abs(d0[0]);
4166 
4167  pointingAngle = selectedMother->CosPointingAngle();
4168  impactProduct = selectedMother->Prodd0d0();
4169  impactProductXY = TMath::Abs(selectedMother->ImpParXY());
4170  invariantMassMother = selectedMother->InvMass(2, prongs);
4171  dcaMother = selectedMother->GetDCA();
4172  vertexDistance = vertexMother->DistanceToVertex(primaryVertex);
4173  d0firstTrack = TMath::Abs(selectedMother->Getd0Prong(0));
4174  d0secondTrack = TMath::Abs(selectedMother->Getd0Prong(1));
4175  normDecayLength = selectedMother->NormalizedDecayLength();
4176 
4177  Double_t pseudoProperDecayLength = ((vertexMother->GetX() - primaryVertex->GetX()) * selectedMother->Px() / TMath::Abs(selectedMother->Pt())) + ((vertexMother->GetY() - primaryVertex->GetY()) * selectedMother->Py() / TMath::Abs(selectedMother->Pt()));
4178  Double_t pseudoProperDecayTime = pseudoProperDecayLength * pdgMassMother / ptMother;
4179  decayTime = vertexDistance / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother / (momentumMother * momentumMother)) + 1)));
4180 
4181  Double_t phi = selectedMother->Phi();
4182  Double_t theta = selectedMother->Theta();
4183  Double_t covMatrix[21];
4184  selectedMother->GetCovarianceXYZPxPyPz(covMatrix);
4185 
4186  Double_t cp = TMath::Cos(phi);
4187  Double_t sp = TMath::Sin(phi);
4188  Double_t ct = TMath::Cos(theta);
4189  Double_t st = TMath::Sin(theta);
4190 
4191  Double_t errorMomentum = covMatrix[9] * cp * cp * ct * ct // GetCovPxPx
4192  + covMatrix[13] * 2.*cp * sp * ct * ct // GetCovPxPy
4193  + covMatrix[18] * 2.*cp * ct * st // GetCovPxPz
4194  + covMatrix[14] * sp * sp * ct * ct // GetCovPyPy
4195  + covMatrix[19] * 2.*sp * ct * st // GetCovPyPz
4196  + covMatrix[20] * st * st; // GetCovPzPz
4197  Double_t normalizedDecayTime = selectedMother->NormalizedDecayLength() / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother * errorMomentum * errorMomentum / (momentumMother * momentumMother)) + 1)));
4198 
4199  Double_t vertexMotherX = vertexMother->GetX();
4200  Double_t vertexMotherY = vertexMother->GetY();
4201  Double_t vertexMotherZ = vertexMother->GetZ();
4202 
4203  Double_t cosPointingAngleXY = selectedMother->CosPointingAngleXY();
4204  Double_t distanceXYToVertex = vertexMother->DistanceXYToVertex(primaryVertex);
4205  Double_t normalizedDecayLengthXY = selectedMother->NormalizedDecayLengthXY();
4206 
4207  ((TH1F*)fMotherHistogramArray[motherType][histType][0])->Fill(ptMother);
4208  ((TH1F*)fMotherHistogramArray[motherType][histType][1])->Fill(ptFirstDaughter);
4209  ((TH1F*)fMotherHistogramArray[motherType][histType][2])->Fill(ptSecondDaughter);
4210  ((TH1F*)fMotherHistogramArray[motherType][histType][3])->Fill(etaMother);
4211  ((TH1F*)fMotherHistogramArray[motherType][histType][4])->Fill(phiMother);
4212  ((TH1F*)fMotherHistogramArray[motherType][histType][5])->Fill(d0Mother);
4213  ((TH1F*)fMotherHistogramArray[motherType][histType][6])->Fill(d0firstTrack);
4214  ((TH1F*)fMotherHistogramArray[motherType][histType][7])->Fill(d0secondTrack);
4215  ((TH1F*)fMotherHistogramArray[motherType][histType][8])->Fill(pointingAngle);
4216  ((TH1F*)fMotherHistogramArray[motherType][histType][9])->Fill(impactProduct);
4217  ((TH1F*)fMotherHistogramArray[motherType][histType][10])->Fill(impactProductXY);
4218  ((TH1F*)fMotherHistogramArray[motherType][histType][11])->Fill(invariantMassMother);
4219  ((TH1F*)fMotherHistogramArray[motherType][histType][12])->Fill(invmassDelta);
4220  ((TH1F*)fMotherHistogramArray[motherType][histType][13])->Fill(dcaMother);
4221  ((TH1F*)fMotherHistogramArray[motherType][histType][14])->Fill(vertexDistance);
4222  ((TH1F*)fMotherHistogramArray[motherType][histType][15])->Fill(normDecayLength);
4223  ((TH1F*)fMotherHistogramArray[motherType][histType][16])->Fill(pseudoProperDecayTime);
4224  ((TH1F*)fMotherHistogramArray[motherType][histType][17])->Fill(decayTime);
4225  ((TH1F*)fMotherHistogramArray[motherType][histType][18])->Fill(normalizedDecayTime);
4226  ((TH1F*)fMotherHistogramArray[motherType][histType][19])->Fill(angleMotherFirstDaughter);
4227  ((TH1F*)fMotherHistogramArray[motherType][histType][20])->Fill(angleMotherSecondDaughter);
4228  ((TH1F*)fMotherHistogramArray[motherType][histType][21])->Fill(angleBetweenBothDaughters);
4229  ((TH1F*)fMotherHistogramArray[motherType][histType][22])->Fill(cosThetaStar);
4230 
4231  ((TH1F*)fMotherHistogramArray[motherType][histType][23])->Fill(vertexMotherX);
4232  ((TH1F*)fMotherHistogramArray[motherType][histType][24])->Fill(vertexMotherY);
4233  ((TH1F*)fMotherHistogramArray[motherType][histType][25])->Fill(vertexMotherZ);
4234 
4235  ((TH1F*)fMotherHistogramArray[motherType][histType][36])->Fill(TMath::Abs(dd0pr1));
4236  ((TH1F*)fMotherHistogramArray[motherType][histType][37])->Fill(TMath::Abs(dd0pr2));
4237  ((TH1F*)fMotherHistogramArray[motherType][histType][38])->Fill(TMath::Abs(dd0max));
4238  ((TH1F*)fMotherHistogramArray[motherType][histType][39])->Fill(TMath::Abs(dd0min));
4239 
4240  ((TH1F*)fMotherHistogramArray[motherType][histType][40])->Fill(cosPointingAngleXY);
4241  ((TH1F*)fMotherHistogramArray[motherType][histType][41])->Fill(distanceXYToVertex);
4242  ((TH1F*)fMotherHistogramArray[motherType][histType][42])->Fill(normalizedDecayLengthXY);
4243  ((TH1F*)fMotherHistogramArray[motherType][histType][43])->Fill(vertexMother->GetChi2());
4244  ((TH1F*)fMotherHistogramArray[motherType][histType][44])->Fill(vertexMother->GetChi2perNDF());
4245 
4246 
4247  //we fill the 2D histograms
4248  Int_t nFirst = 0;
4249  Int_t nSecond = 1;
4250  Int_t nVariables = 10;
4251  Int_t nHistograms = nVariables * (nVariables - 1) / 2;
4252  for (Int_t k = 0; k < nHistograms; ++k)
4253  {
4254  Double_t firstVariable = 0.0;
4255  Double_t secondVariable = 0.0;
4256 
4257  if(nFirst==0) firstVariable = d0firstTrack;
4258  if(nFirst==1) firstVariable = d0secondTrack;
4259  if(nFirst==2) firstVariable = d0Mother;
4260  if(nFirst==3) firstVariable = pointingAngle;
4261  if(nFirst==4) firstVariable = impactProduct;
4262  if(nFirst==5) firstVariable = impactProductXY;
4263  if(nFirst==6) firstVariable = vertexDistance;
4264  if(nFirst==7) firstVariable = normDecayLength;
4265  if(nFirst==8) firstVariable = cosPointingAngleXY;
4266  if(nFirst==9) firstVariable = distanceXYToVertex;
4267  if(nFirst==10) firstVariable = normalizedDecayLengthXY;
4268 
4269  if(nSecond==0) secondVariable = d0firstTrack;
4270  if(nSecond==1) secondVariable = d0secondTrack;
4271  if(nSecond==2) secondVariable = d0Mother;
4272  if(nSecond==3) secondVariable = pointingAngle;
4273  if(nSecond==4) secondVariable = impactProduct;
4274  if(nSecond==5) secondVariable = impactProductXY;
4275  if(nSecond==6) secondVariable = vertexDistance;
4276  if(nSecond==7) secondVariable = normDecayLength;
4277  if(nSecond==8) secondVariable = cosPointingAngleXY;
4278  if(nSecond==9) secondVariable = distanceXYToVertex;
4279  if(nSecond==10) secondVariable = normalizedDecayLengthXY;
4280 
4281  ((TH2F*)fMotherHistogramArray2D[motherType][histType][k])->Fill(firstVariable,secondVariable);
4282 
4283  nSecond++;
4284  if(nSecond>nVariables)
4285  {
4286  nFirst++;
4287  nSecond = nFirst + 1;
4288  }
4289  }
4290 
4291 
4292  if (motherType == 1) {
4293  motherType = motherType - 1;
4294 
4295  AliAODRecoDecay* trackD0 = (AliAODRecoDecay*)selectedMother->GetDaughter(1);
4296  AliAODTrack * pionD0 = (AliAODTrack*)trackD0->GetDaughter(0);
4297  AliAODTrack * kaonD0 = (AliAODTrack*)trackD0->GetDaughter(1);
4298 
4299  AliAODVertex * vertexBPlus = vertexMother;
4300  AliAODVertex * vertexD0 = trackD0->GetSecondaryVtx();
4301  vertexDistance = TMath::Abs(vertexBPlus->DistanceToVertex(vertexD0));
4302  pdgMassMother = TDatabasePDG::Instance()->GetParticle(421)->Mass();
4303 
4304  AliExternalTrackParam pionD0Track;
4305  AliExternalTrackParam kaonD0Track;
4306 
4307  Double_t d0z0[2], covd0z0[3], d0[2];
4308 
4309  pionD0Track.CopyFromVTrack(pionD0);
4310  pionD0Track.PropagateToDCA(vertexBPlus, bz, 100., d0z0, covd0z0);
4311  d0[0] = d0z0[0];
4312 
4313  kaonD0Track.CopyFromVTrack(kaonD0);
4314  kaonD0Track.PropagateToDCA(vertexBPlus, bz, 100., d0z0, covd0z0);
4315  d0[1] = d0z0[0];
4316 
4317  AliExternalTrackParam D0Track;
4318  D0Track.CopyFromVTrack(trackD0);
4319  Double_t d0z0D0[2], covd0z0D0[3], d0D0;
4320  D0Track.PropagateToDCA(vertexBPlus, bz, 100., d0z0D0, covd0z0D0);
4321  d0D0 = d0z0D0[0];
4322 
4323  Double_t impactProductToBPlus = d0[0] * d0[1];
4324  Double_t impactProductXYToBPlus = trackD0->ImpParXY(vertexBPlus);
4325 
4326  Double_t momentumMother = trackD0->P();
4327  Double_t pointingAngleToBPlus = trackD0->CosPointingAngle(vertexBPlus);
4328  Double_t d0FirstDaughterToBPlus = TMath::Abs(d0[0]);
4329  Double_t d0trackD0ToBPlus = TMath::Abs(d0[1]);
4330  Double_t normDecayLengthToBPlus = trackD0->NormalizedDecayLength(vertexBPlus);
4331 
4332  Double_t pseudoProperDecayLength = ((vertexD0->GetX() - vertexBPlus->GetX()) * trackD0->Px() / TMath::Abs(trackD0->Pt())) + ((vertexD0->GetY() - vertexBPlus->GetY()) * trackD0->Py() / TMath::Abs(trackD0->Pt()));
4333  Double_t pseudoProperDecayTimeToBPlus = pseudoProperDecayLength * pdgMassMother / ptMother;
4334  Double_t DecayTimeToBPlus = vertexDistance / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother / (momentumMother * momentumMother)) + 1)));
4335 
4336  Double_t phi = trackD0->Phi();
4337  Double_t theta = trackD0->Theta();
4338  Double_t covMatrix[21];
4339  trackD0->GetCovarianceXYZPxPyPz(covMatrix);
4340 
4341  Double_t cp = TMath::Cos(phi);
4342  Double_t sp = TMath::Sin(phi);
4343  Double_t ct = TMath::Cos(theta);
4344  Double_t st = TMath::Sin(theta);
4345 
4346  Double_t errorMomentum = covMatrix[9] * cp * cp * ct * ct // GetCovPxPx
4347  + covMatrix[13] * 2.*cp * sp * ct * ct // GetCovPxPy
4348  + covMatrix[18] * 2.*cp * ct * st // GetCovPxPz
4349  + covMatrix[14] * sp * sp * ct * ct // GetCovPyPy
4350  + covMatrix[19] * 2.*sp * ct * st // GetCovPyPz
4351  + covMatrix[20] * st * st; // GetCovPzPz
4352  Double_t normDecayTimeToBPlus = trackD0->NormalizedDecayLength(vertexBPlus) / (299792458 * TMath::Sqrt(1 / ((pdgMassMother * pdgMassMother * errorMomentum * errorMomentum / (momentumMother * momentumMother)) + 1)));
4353 
4354  ((TH1F*)fMotherHistogramArray[motherType][histType][26])->Fill(pointingAngleToBPlus);
4355  ((TH1F*)fMotherHistogramArray[motherType][histType][27])->Fill(d0D0);
4356  ((TH1F*)fMotherHistogramArray[motherType][histType][28])->Fill(d0FirstDaughterToBPlus);
4357  ((TH1F*)fMotherHistogramArray[motherType][histType][29])->Fill(d0trackD0ToBPlus);
4358  ((TH1F*)fMotherHistogramArray[motherType][histType][30])->Fill(impactProductToBPlus);
4359  ((TH1F*)fMotherHistogramArray[motherType][histType][31])->Fill(impactProductXYToBPlus);
4360  ((TH1F*)fMotherHistogramArray[motherType][histType][32])->Fill(normDecayLengthToBPlus);
4361  ((TH1F*)fMotherHistogramArray[motherType][histType][33])->Fill(pseudoProperDecayTimeToBPlus);
4362  ((TH1F*)fMotherHistogramArray[motherType][histType][34])->Fill(DecayTimeToBPlus);
4363  ((TH1F*)fMotherHistogramArray[motherType][histType][35])->Fill(normDecayTimeToBPlus);
4364 
4365 
4366 
4367  }
4368  return;
4369 }
4370 //-------------------------------------------------------------------------------------
4371 void AliAnalysisTaskSEBPlustoD0Pi::TwoTrackCombinationInfo(AliExternalTrackParam * firstTrack, AliExternalTrackParam * secondTrack, AliAODVertex * primaryVertex, Double_t bz, Bool_t isDesiredCandidate, TString histogram_name, UInt_t prongs[2]) {
4372 
4373  // we calculate the vertex position
4374  TObjArray daughterTracks;
4375 
4376  daughterTracks.Add(firstTrack);
4377  daughterTracks.Add(secondTrack);
4378 
4379  Double_t dispersion = 0;
4380  AliAODVertex *vertex = RecalculateVertex(primaryVertex, &daughterTracks, bz, dispersion);
4381  if (!vertex) {delete vertex; vertex = NULL; return;}
4382 
4383  Double_t xdummy = 0., ydummy = 0., dca, e[2];
4384  Double_t d0z0[2], covd0z0[3], d0[2], d0err[2];
4385 
4386  firstTrack->PropagateToDCA(vertex, bz, 100., d0z0, covd0z0);
4387  secondTrack->PropagateToDCA(vertex, bz, 100., d0z0, covd0z0);
4388  dca = secondTrack->GetDCA(firstTrack, bz, xdummy, ydummy);
4389 
4390  Double_t px[2], py[2], pz[2];
4391  px[0] = firstTrack->Px();
4392  py[0] = firstTrack->Py();
4393  pz[0] = firstTrack->Pz();
4394  px[1] = secondTrack->Px();
4395  py[1] = secondTrack->Py();
4396  pz[1] = secondTrack->Pz();
4397 
4398  firstTrack->PropagateToDCA(primaryVertex, bz, 100., d0z0, covd0z0);
4399  d0[0] = d0z0[0];
4400  d0err[0] = TMath::Sqrt(covd0z0[0]);
4401  secondTrack->PropagateToDCA(primaryVertex, bz, 100., d0z0, covd0z0);
4402  d0[1] = d0z0[0];
4403  d0err[1] = TMath::Sqrt(covd0z0[0]);
4404 
4405  AliAODRecoDecayHF2Prong * track = new AliAODRecoDecayHF2Prong(vertex, px, py, pz, d0, d0err, dca);
4406  if (!track)
4407  {
4408  delete vertex; vertex = NULL;
4409  delete track; track = NULL;
4410  return;
4411  }
4412 
4413  Double_t invariantMass = track->InvMass(2, prongs);
4414 
4415  TString fill = histogram_name;
4416  if (isDesiredCandidate) fill += "_MC";
4417  ((TH1F*)(fOutputBPlusMC->FindObject(fill)))->Fill(invariantMass);
4418 
4419  delete vertex; vertex = NULL;
4420  delete track; track = NULL;
4421  return;
4422 }
4423 //-------------------------------------------------------------------------------------
4424 void AliAnalysisTaskSEBPlustoD0Pi::ThreeTrackCombinationInfo(AliExternalTrackParam * firstTrack, AliExternalTrackParam * secondTrack, AliExternalTrackParam * thirdTrack, AliAODVertex * primaryVertex, Double_t bz, Bool_t isDesiredCandidate, TString histogram_name, UInt_t prongs[3]) {
4425 
4426  // we calculate the vertex position
4427  TObjArray daughterTracks;
4428 
4429  daughterTracks.Add(firstTrack);
4430  daughterTracks.Add(secondTrack);
4431  daughterTracks.Add(thirdTrack);
4432 
4433  Double_t dispersion = 0;
4434  AliAODVertex *vertex = RecalculateVertex(primaryVertex, &daughterTracks, bz, dispersion);
4435  if (!vertex) {delete vertex; vertex = NULL; return;}
4436 
4437  Double_t xdummy = 0., ydummy = 0., dca[3];
4438  Double_t d0z0[2], covd0z0[3], d0[3], d0err[3];
4439 
4440  firstTrack->PropagateToDCA(vertex, bz, 100., d0z0, covd0z0);
4441  secondTrack->PropagateToDCA(vertex, bz, 100., d0z0, covd0z0);
4442  thirdTrack->PropagateToDCA(vertex, bz, 100., d0z0, covd0z0);
4443 
4444  dca[0] = firstTrack->GetDCA(secondTrack, bz, xdummy, ydummy);
4445  dca[1] = firstTrack->GetDCA(thirdTrack, bz, xdummy, ydummy);
4446  dca[2] = secondTrack->GetDCA(thirdTrack, bz, xdummy, ydummy);
4447 
4448  Double_t px[3], py[3], pz[3];
4449  px[0] = firstTrack->Px();
4450  py[0] = firstTrack->Py();
4451  pz[0] = firstTrack->Pz();
4452  px[1] = secondTrack->Px();
4453  py[1] = secondTrack->Py();
4454  pz[1] = secondTrack->Pz();
4455  px[2] = thirdTrack->Px();
4456  py[2] = thirdTrack->Py();
4457  pz[2] = thirdTrack->Pz();
4458 
4459  firstTrack->PropagateToDCA(primaryVertex, bz, 100., d0z0, covd0z0);
4460  d0[0] = d0z0[0];
4461  d0err[0] = TMath::Sqrt(covd0z0[0]);
4462  secondTrack->PropagateToDCA(primaryVertex, bz, 100., d0z0, covd0z0);
4463  d0[1] = d0z0[0];
4464  d0err[1] = TMath::Sqrt(covd0z0[0]);
4465  thirdTrack->PropagateToDCA(primaryVertex, bz, 100., d0z0, covd0z0);
4466  d0[2] = d0z0[0];
4467  d0err[2] = TMath::Sqrt(covd0z0[0]);
4468 
4469  // Double_t pos[3]; primaryVertex->GetXYZ(pos);
4470  // Double_t dist12=TMath::Sqrt((vertex12->GetX()-pos[0])*(vertex12->GetX()-pos[0])+(vertex12->GetY()-pos[1])*(vertex12->GetY()-pos[1])+(vertex12->GetZ()-pos[2])*(vertex12->GetZ()-pos[2]));
4471  // Double_t dist23=TMath::Sqrt((vertex23->GetX()-pos[0])*(vertex23->GetX()-pos[0])+(vertex23->GetY()-pos[1])*(vertex23->GetY()-pos[1])+(vertex23->GetZ()-pos[2])*(vertex23->GetZ()-pos[2]));
4472  Short_t charge = (Short_t)(firstTrack->Charge() + secondTrack->Charge() + thirdTrack->Charge());
4473 
4474  AliAODRecoDecayHF3Prong * track = new AliAODRecoDecayHF3Prong(vertex, px, py, pz, d0, d0err, dca, dispersion, 0.0, 0.0, charge); //dist12,dist23,charge);
4475  if (!track)
4476  {
4477  delete vertex; vertex = NULL;
4478  delete track; track = NULL;
4479  return;
4480  }
4481 
4482  Double_t invariantMass = track->InvMass(3, prongs);
4483 
4484  TString fill = histogram_name;
4485  if (isDesiredCandidate) fill += "_MC";
4486  ((TH1F*)(fOutputBPlusMC->FindObject(fill)))->Fill(invariantMass);
4487 
4488  delete vertex; vertex = NULL;
4489  delete track; track = NULL;
4490  return;
4491 }
4492 //-------------------------------------------------------------------------------------
4493 Int_t AliAnalysisTaskSEBPlustoD0Pi::MatchCandidateToMonteCarlo(Int_t pdgabs, AliAODRecoDecayHF2Prong * candidate, TClonesArray *mcArray, TMatrix * BPlustoD0PiLabelMatrix) const
4494 {
4495  //
4496  // Check if this candidate is matched to a MC signal
4497  // If no, return -1
4498  // If yes, return label (>=0) of the AliAODMCParticle
4499 
4500 
4501  // Check number of daughters
4502  Int_t ndg = candidate->GetNDaughters();
4503  if(!ndg) { AliError("No daughters available"); return -1;}
4504  if(ndg != 2) return -1;
4505 
4506  // loop on daughters and write the labels
4507  Int_t dgLabels[2] = {-1};
4508  Int_t pdgDg[2] = {0};
4509  Int_t signalPosition = -1;
4510  if(pdgabs==421)
4511  {
4512  AliAODTrack *trk0 = (AliAODTrack*)candidate->GetDaughter(0);
4513  dgLabels[0] = trk0->GetLabel();
4514  AliAODTrack *trk1 = (AliAODTrack*)candidate->GetDaughter(1);
4515  dgLabels[1] = trk1->GetLabel();
4516  pdgDg[0] = 211; pdgDg[1] = 321;
4517  signalPosition = 3;
4518  }
4519  else if(pdgabs==521)
4520  {
4521  AliAODTrack *trk0 = (AliAODTrack*)candidate->GetDaughter(0);
4522  dgLabels[0] = trk0->GetLabel();
4523  dgLabels[1] = MatchCandidateToMonteCarlo(421,(AliAODRecoDecayHF2Prong*)candidate->GetDaughter(1), mcArray, BPlustoD0PiLabelMatrix);
4524  pdgDg[0] = 211; pdgDg[1] = 421;
4525  signalPosition = 4;
4526  }
4527  else
4528  {
4529  std::cout << "Wrong pdg supplied for function to match candidate to monte carlo signal." << std::endl;
4530  return -1;
4531  }
4532  if(dgLabels[0]==-1) return -1;
4533  if(dgLabels[1]==-1) return -1;
4534 
4535 
4536  Int_t labMom[2]={0,0};
4537  Int_t i,j,lab,labMother,pdgMother,pdgPart;
4538  AliAODMCParticle *part=0;
4539  AliAODMCParticle *mother=0;
4540  Double_t pxSumDgs=0.,pySumDgs=0.,pzSumDgs=0.;
4541  Bool_t pdgUsed[2]={kFALSE,kFALSE};
4542 
4543  // loop on daughter labels
4544  for(i=0; i<ndg; i++)
4545  {
4546  labMom[i]=-1;
4547  lab = TMath::Abs(dgLabels[i]);
4548  if(lab<0)
4549  {
4550  printf("daughter with negative label %d\n",lab);
4551  return -1;
4552  }
4553  part = (AliAODMCParticle*)mcArray->At(lab);
4554  if(!part)
4555  {
4556  printf("no MC particle\n");
4557  return -1;
4558  }
4559 
4560  // check the PDG of the daughter
4561  pdgPart=TMath::Abs(part->GetPdgCode());
4562  for(j=0; j<ndg; j++)
4563  {
4564  if(!pdgUsed[j] && pdgPart==pdgDg[j])
4565  {
4566  pdgUsed[j]=kTRUE;
4567  break;
4568  }
4569  }
4570 
4571 
4572  mother = part;
4573  while(mother->GetMother()>=0)
4574  {
4575  labMother=mother->GetMother();
4576  mother = (AliAODMCParticle*)mcArray->At(labMother);
4577  if(!mother)
4578  {
4579  printf("no MC mother particle\n");
4580  break;
4581  }
4582  pdgMother = TMath::Abs(mother->GetPdgCode());
4583  if(pdgMother==pdgabs)
4584  {
4585  labMom[i]=labMother;
4586  // keep sum of daughters' momenta, to check for mom conservation
4587  pxSumDgs += part->Px();
4588  pySumDgs += part->Py();
4589  pzSumDgs += part->Pz();
4590  break;
4591  }
4592  else break;
4593  }
4594  if(labMom[i]==-1) return -1; // mother PDG not ok for this daughter
4595  } // end loop on daughters
4596 
4597  // check if the candidate is signal
4598  labMother=labMom[0];
4599  // all labels have to be the same and !=-1
4600  for(i=0; i<ndg; i++)
4601  {
4602  if(labMom[i]==-1) return -1;
4603  if(labMom[i]!=labMother) return -1;
4604  }
4605 
4606  // check that all daughter PDGs are matched
4607  for(i=0; i<ndg; i++)
4608  {
4609  if(pdgUsed[i]==kFALSE) return -1;
4610  }
4611 
4612  // Check for mom conservation
4613  mother = (AliAODMCParticle*)mcArray->At(labMother);
4614  Double_t pxMother = mother->Px();
4615  Double_t pyMother = mother->Py();
4616  Double_t pzMother = mother->Pz();
4617  // within 0.5%
4618  if((TMath::Abs(pxMother-pxSumDgs)/(TMath::Abs(pxMother)+1.e-13)) > 0.05 &&
4619  (TMath::Abs(pyMother-pySumDgs)/(TMath::Abs(pyMother)+1.e-13)) > 0.05 &&
4620  (TMath::Abs(pzMother-pzSumDgs)/(TMath::Abs(pzMother)+1.e-13)) > 0.05)
4621  {
4622  return -1;
4623  }
4624 
4625  // Check if label matches a signal label
4626  // Int_t bIsSignal = kFALSE;
4627  // TMatrix &particleMatrix = *BPlustoD0PiLabelMatrix;
4628  // for (Int_t k = 0; k < BPlustoD0PiLabelMatrix->GetNrows(); ++k)
4629  // {
4630  // if(labMother == (Int_t)particleMatrix(k,signalPosition))
4631  // {
4632  // bIsSignal = kTRUE;
4633  // break;
4634  // }
4635  // }
4636  // if(!bIsSignal) return -1;
4637 
4638  return labMother;
4639 }
4640 //-------------------------------------------------------------------------------------
4641 Int_t AliAnalysisTaskSEBPlustoD0Pi::IsTrackInjected(AliAODTrack *part,AliAODMCHeader *header,TClonesArray *arrayMC){
4642 
4644 
4645  Int_t lab=part->GetLabel();
4646  if(lab<0) {delete ggg; ggg = nullptr; return 1;} //
4647  TString nameGen = ggg->GetGenerator(lab,header);
4648  TString empty="";
4649  Int_t countControl =0;
4650  while(nameGen.IsWhitespace()){
4651  AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);
4652  if(!mcpart){
4653  printf("AliVertexingHFUtils::IsTrackInjected - BREAK: No valid AliAODMCParticle at label %i\n",lab);
4654  break;
4655  }
4656  Int_t mother = mcpart->GetMother();
4657  if(mother<0){
4658  // printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Reached primary particle without valid mother\n");
4659  break;
4660  }
4661  lab=mother;
4662  nameGen = ggg->GetGenerator(mother,header);
4663  countControl++;
4664  if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops
4665  printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Protection from infinite loop active\n");
4666  break;
4667  }
4668  }
4669  if(nameGen.IsWhitespace() || nameGen.Contains("ijing")){delete ggg; ggg = nullptr; return 0;}
4670 
4671  delete ggg; ggg = nullptr;
4672  return 1;
4673 }
4674 //-------------------------------------------------------------------------------------
4675 Bool_t AliAnalysisTaskSEBPlustoD0Pi::IsCandidateInjected(AliAODRecoDecayHF2Prong *selectedBPlus, AliAODMCHeader *header,TClonesArray *arrayMC){
4676 
4677  AliAODTrack* selectedBPlusPion = (AliAODTrack*)selectedBPlus->GetDaughter(0);
4678  AliAODRecoDecayHF2Prong* selectedD0 = (AliAODRecoDecayHF2Prong*)selectedBPlus->GetDaughter(1);
4679 
4680  AliAODTrack* selectedD0FirstDaughter = (AliAODTrack*)selectedD0->GetDaughter(0);
4681  AliAODTrack* selectedD0SecondDaughter = (AliAODTrack*)selectedD0->GetDaughter(1);
4682 
4683  if(IsTrackInjected(selectedBPlusPion,header,arrayMC)) return kTRUE;
4684  if(IsTrackInjected(selectedD0FirstDaughter,header,arrayMC)) return kTRUE;
4685  if(IsTrackInjected(selectedD0SecondDaughter,header,arrayMC)) return kTRUE;
4686 
4687  return kFALSE;
4688 }
Int_t charge
Double_t NormalizedDecayLengthXY() const
Double_t NormalizedDecayLength() const
Int_t IsTrackInjected(AliAODTrack *part, AliAODMCHeader *header, TClonesArray *arrayMC)
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel, AliAODEvent *aod, Bool_t *bCutArray)
double Double_t
Definition: External.C:58
ULong64_t GetTriggerMask()
Definition: AliRDHFCuts.h:69
Double_t DeltaInvMassBPlusKpipi(AliAODRecoDecayHF2Prong *BPlus) const
Float_t * GetPtBinLimitsD0forD0ptbin() const
Definition: External.C:236
void Getd0MeasMinusExpProng(Int_t ip, Double_t magf, Double_t &d0diff, Double_t &errd0diff) const
Int_t GetnSigmaTOF(AliAODTrack *track, Int_t species, Double_t &sigma) const
void GetSoftSelectionArrayITSD0SecondDaughter(Bool_t array[7]=0)
Bool_t D0FirstDaughterSelection(AliAODTrack *aodTrack, AliAODVertex *primaryVertex, Double_t bz, TClonesArray *mcTrackArray, TMatrix *B0toDStarPiLabelMatrix, AliAODMCHeader *header)
Int_t GetnSigmaTPC(AliAODTrack *track, Int_t species, Double_t &sigma) const
Double_t bz
Double_t ImpParXY() const
char Char_t
Definition: External.C:18
static TString GetGenerator(Int_t label, AliAODMCHeader *header)
Int_t GetWhyRejection() const
Definition: AliRDHFCuts.h:315
Double_t CosPointingAngleXY() const
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
void BPlustoD0PiSignalTracksInMC(TClonesArray *mcTrackArray, AliAODEvent *aodevent, TMatrix *BPlustoD0PiLabelMatrix, TList *listout)
Int_t IsD0forD0ptbinSelected(TObject *obj, Int_t selectionLevel, AliAODEvent *aod, Bool_t *bCutArray)
AliAODPidHF * GetPidHF() const
Definition: AliRDHFCuts.h:247
virtual void Terminate(Option_t *option)
AliAODVertex * RecalculateVertex(const AliVVertex *primary, TObjArray *tracks, Double_t bField, Double_t dispersion)
void BPlusPionSelection(AliAODEvent *aodEvent, AliAODVertex *primaryVertex, Double_t bz, TClonesArray *mcTrackArray, TMatrix *BPlustoD0PiLabelMatrix, AliAODMCHeader *header)
Bool_t D0SecondDaughterSelection(AliAODTrack *aodTrack, AliAODVertex *primaryVertex, Double_t bz, TClonesArray *mcTrackArray, TMatrix *B0toDStarPiLabelMatrix, AliAODMCHeader *header)
void TwoTrackCombinationInfo(AliExternalTrackParam *firstTrack, AliExternalTrackParam *secondTrack, AliAODVertex *primaryVertex, Double_t bz, Bool_t isDesiredCandidate, TString histogram_name, UInt_t prongs[2])
Int_t PtBinD0forD0ptbin(Double_t pt) const
int Int_t
Definition: External.C:63
Definition: External.C:204
unsigned int UInt_t
Definition: External.C:33
void GetHardSelectionArrayITSD0FirstDaughter(Bool_t array[7]=0)
float Float_t
Definition: External.C:68
virtual Int_t SelectPID(AliAODTrack *track, Int_t type)
Int_t MatchCandidateToMonteCarlo(Int_t pdgabs, AliAODRecoDecayHF2Prong *candidate, TClonesArray *mcArray, TMatrix *B0toDStarPiLabelMatrix) const
void SetProngIDs(Int_t nIDs, UShort_t *id)
void GetHardSelectionArrayITSD0SecondDaughter(Bool_t array[7]=0)
void ThreeTrackCombinationInfo(AliExternalTrackParam *firstTrack, AliExternalTrackParam *secondTrack, AliExternalTrackParam *thirdTrack, AliAODVertex *primaryVertex, Double_t bz, Bool_t isDesiredCandidate, TString histogram_name, UInt_t prongs[3])
void SetPrimaryVtxRef(TObject *vtx)
primary vertex
short Short_t
Definition: External.C:23
void GetSoftSelectionArrayITSD0FirstDaughter(Bool_t array[7]=0)
virtual void UserCreateOutputObjects()
Implementation of interface methods.
Bool_t IsEventSelected(AliVEvent *event)
void FillBPlusHistograms(AliAODRecoDecayHF2Prong *selectedMother, AliAODVertex *primaryVertex, Double_t bz, Int_t motherType, Int_t histType)
void FillD0Histograms(AliAODRecoDecayHF2Prong *selectedMother, AliAODVertex *primaryVertex, Double_t bz, Int_t motherType, Int_t histType, Int_t pdgCodeMother=-1)
void FillFinalTrackHistograms(AliAODRecoDecayHF2Prong *selectedBPlus, AliAODVertex *primaryVertex, Double_t bz, Bool_t isDesiredCandidate, TClonesArray *mcTrackArray)
void GetSoftSelectionArrayITSBPlusPion(Bool_t array[7]=0)
Float_t * GetPtBinLimits() const
Definition: AliRDHFCuts.h:248
Bool_t IsCandidateInjected(AliAODRecoDecayHF2Prong *part, AliAODMCHeader *header, TClonesArray *arrayMC)
void GetHardSelectionArrayITSBPlusPion(Bool_t array[7]=0)
unsigned short UShort_t
Definition: External.C:28
const char Option_t
Definition: External.C:48
TH1F * fDaughterHistogramArray[4][6][15]
[fnPtBinsD0forD0ptbinLimits]
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:249
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
virtual void UserExec(Option_t *option)
Int_t GetNPtBinsD0forD0ptbin() const
Int_t PtBin(Double_t pt) const
void BPlusSelection(AliAODEvent *aodEvent, AliAODVertex *primaryVertex, Double_t bz, TClonesArray *mcTrackArray, TMatrix *BPlustoD0PiLabelMatrix, TClonesArray *D0TracksFromFriendFile, AliAODMCHeader *header)
void D0Selection(AliAODEvent *aodEvent, AliAODVertex *primaryVertex, Double_t bz, TClonesArray *mcTrackArray, TMatrix *BPlustoD0PiLabelMatrix, TClonesArray *D0TracksFromFriendFile, AliAODMCHeader *header)
Class with functions useful for different D2H analyses //.