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