AliPhysics  5a28df1 (5a28df1)
AliAnalysisTaskSEB0toDStarPi.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 (B0 -> DStar pi -> D0 pi pi -> K pi pi pi) Analysis
21 //
22 //
23 // Cuts are centralized in AliRDHFCutsB0toDStarPi
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 // Based on the DStartoKpipi macro by:
32 //
33 // Author A.Grelli
34 // ERC-QGP Utrecht University - a.grelli@uu.nl,
35 // Author Y.Wang
36 // University of Heidelberg - yifei@physi.uni-heidelberg.de
37 // Author C.Ivan
38 // ERC-QGP Utrecht University - c.ivan@uu.nl,
39 //
40 //-----------------------------------------------------------------------
41 
42 #include <TSystem.h>
43 #include <TParticle.h>
44 #include <TH1I.h>
45 #include "TROOT.h"
46 #include <TDatabasePDG.h>
47 #include <AliAnalysisDataSlot.h>
48 #include <AliAnalysisDataContainer.h>
49 #include "AliRDHFCutsB0toDStarPi.h"
50 #include "AliStack.h"
51 #include "AliMCEvent.h"
52 #include "AliAnalysisManager.h"
53 #include "AliAODMCHeader.h"
54 #include "AliAODHandler.h"
55 #include "AliLog.h"
56 #include "AliVertex.h"
57 #include "AliVVertex.h"
58 #include "AliESDVertex.h"
59 #include "AliAODVertex.h"
60 #include "AliVertexerTracks.h"
61 #include "AliExternalTrackParam.h"
62 #include "AliNeutralTrackParam.h"
63 #include "AliAODRecoDecay.h"
64 #include "AliAODRecoDecayHF.h"
67 #include "AliAnalysisVertexingHF.h"
68 #include "AliESDtrack.h"
69 #include "AliAODMCParticle.h"
71 #include "AliAODEvent.h"
73 #include "AliAODInputHandler.h"
74 #include <vector>
75 #include <TMatrix.h>
76 #include <TVector3.h>
77 #include <TArrayI.h>
78 #include <bitset>
79 
80 // #include "TObjectTable.h"
81 
85 
86 //__________________________________________________________________________
89  fEvents(0),
90  fUseMCInfo(kFALSE),
91  fOutput(0),
92  fOutputD0Pion(0),
93  fOutputD0Kaon(0),
94  fOutputDStarPion(0),
95  fOutputB0Pion(0),
96  fOutputD0(0),
97  fOutputDStar(0),
98  fOutputB0(0),
99  fOutputD0_D0Pt(0),
100  fOutputD0_DStarPt(0),
101  fOutputDStar_DStarPt(0),
102  fOutputB0MC(0),
103  fCuts(0),
104  fQuickSignalAnalysis(0),
105  fGetCutInfo(0),
106  fCEvents(0),
107  fCounter(0),
108  fDStarPionTracks(0),
109  fB0PionTracks(0),
110  fD0Tracks(0),
111  fShowMask(0),
112  fShowRejection(0),
113  fnPtBins(0),
114  fnPtBinsD0forD0ptbin(0),
115  fnPtBinsD0forDStarptbin(0),
116  fnPtBinsDStarforDStarptbin(0),
117  fnPtBinLimits(0),
118  fnPtBinsD0forD0ptbinLimits(0),
119  fnPtBinsD0forDStarptbinLimits(0),
120  fnPtBinsDStarforDStarptbinLimits(0),
121  fPtBinLimits(0),
122  fPtBinLimitsD0forD0ptbin(0),
123  fPtBinLimitsD0forDStarptbin(0),
124  fPtBinLimitsDStarforDStarptbin(0),
125  fDaughterHistogramArray(),
126  fDaughterHistogramArray2D(),
127  fMotherHistogramArray(),
128  fMotherHistogramArray2D(),
129  fMotherHistogramArrayExtra()
130 {
131  //
133  //
134 }
135 //___________________________________________________________________________
137  AliAnalysisTaskSE(name),
138  fEvents(0),
139  fUseMCInfo(kFALSE),
140  fOutput(0),
141  fOutputD0Pion(0),
142  fOutputD0Kaon(0),
143  fOutputDStarPion(0),
144  fOutputB0Pion(0),
145  fOutputD0(0),
146  fOutputDStar(0),
147  fOutputB0(0),
148  fOutputD0_D0Pt(0),
151  fOutputB0MC(0),
152  fCuts(0),
154  fGetCutInfo(0),
155  fCEvents(0),
156  fCounter(0),
157  fDStarPionTracks(0),
158  fB0PionTracks(0),
159  fD0Tracks(0),
160  fShowMask(0),
161  fShowRejection(0),
162  fnPtBins(0),
166  fnPtBinLimits(0),
170  fPtBinLimits(0),
179 {
180  //
182  //
183  Info("AliAnalysisTaskSEB0toDStarPi","Calling Constructor");
184 
185  // we prepare vectors and arrays that will save the candidates during the reconstruction
186  fDStarPionTracks = new std::vector<Int_t>;
187  fB0PionTracks = new std::vector<Int_t>;
188  fD0Tracks = new std::vector<Int_t>;
189 
190  // we get the cut file
191  fCuts=cuts;
192 
193  // we get information on the pt bins
198 
199  fnPtBinLimits = fnPtBins + 1;
203 
208 
209  // we create an array of pointers for the histograms. This method is more CPU efficient than looking up each histogram by name.
210  // Automatic option is not yet complete, the arrays are now set manualy in the header file
211 
212  // const Int_t numberOfDaughters = 4;
213  // const Int_t numberOfDaughterHistogramSets = 5;
214  // const Int_t numberOfDaughterHistograms = 15;
215  // const Int_t numberOfDaughterHistograms2D = 6;
216 
217  // Int_t maxHistogramSets = 6 + 2*fnPtBins;
218  // if(2*fnPtBinsD0forD0ptbin > maxHistogramSets) maxHistogramSets = 2*fnPtBinsD0forD0ptbin;
219  // if(2*fnPtBinsD0forDStarptbin > maxHistogramSets) maxHistogramSets = 2*fnPtBinsD0forDStarptbin;
220  // if(2*fnPtBinsDStarforDStarptbin > maxHistogramSets) maxHistogramSets = 2*fnPtBinsDStarforDStarptbin;
221 
222  // const Int_t numberOfOutputs = 6;
223  // const Int_t numberOfMotherHistogramSets = maxHistogramSets;
224  // const Int_t numberOfMotherHistograms = 46;
225  // const Int_t numberOfMotherHistograms2D = 7;
226 
227  // fDaughterHistogramArray = new Int_t*[numberOfDaughters][numberOfDaughterHistogramSets][numberOfDaughterHistograms];
228  // fDaughterHistogramArray2D = new Int_t*[numberOfDaughters][numberOfDaughterHistograms2D];
229  // fMotherHistogramArray = new Int_t*[numberOfOutputs][numberOfMotherHistogramSets][numberOfMotherHistograms];
230  // fMotherHistogramArray2D = new Int_t*[numberOfOutputs][numberOfMotherHistograms2D];
231 
232 
233  DefineOutput(1,TList::Class()); //counters
234  DefineOutput(2,TList::Class()); //All Entries output
235  DefineOutput(3,TList::Class()); //3sigma PID output
236  DefineOutput(4,AliRDHFCutsB0toDStarPi::Class()); //My private output
237  DefineOutput(5,AliNormalizationCounter::Class()); // normalization
238  DefineOutput(6,TList::Class()); // D0 pion output
239  DefineOutput(7,TList::Class()); // D0 kaon output
240  DefineOutput(8,TList::Class()); // DStar pion output
241  DefineOutput(9,TList::Class()); // B0 pion output
242  DefineOutput(10,TList::Class()); // D0 output
243  DefineOutput(11,TList::Class()); // DStar output
244  DefineOutput(12,TList::Class()); // B0 output
245  DefineOutput(13,TList::Class()); // B0 output
246  DefineOutput(14,TList::Class()); // B0 output
247  DefineOutput(15,TList::Class()); // B0 output
248  DefineOutput(16,TList::Class()); // B0MC output
249 }
250 
251 //___________________________________________________________________________
253  //
255  //
256  Info("~AliAnalysisTaskSEB0toDStarPi","Calling Destructor");
257 
258  delete fOutput;
259  delete fOutputD0Pion;
260  delete fOutputD0Kaon;
261  delete fOutputDStarPion;
262  delete fOutputB0Pion;
263  delete fOutputD0;
264  delete fOutputDStar;
265  delete fOutputB0;
266  delete fOutputD0_D0Pt;
267  delete fOutputD0_DStarPt;
268  delete fOutputDStar_DStarPt;
269  delete fOutputB0MC;
270  delete fCuts;
271  delete fCEvents;
272  delete fD0Tracks;
273  delete fDStarPionTracks;
274  delete fB0PionTracks;
275 }
276 //_________________________________________________
278  //
280  //
281 
282  if(fDebug > 1) printf("AliAnalysisTaskSEB0toDStarPi::Init() \n");
283  AliRDHFCutsB0toDStarPi* copyfCuts=new AliRDHFCutsB0toDStarPi(*(static_cast<AliRDHFCutsB0toDStarPi*>(fCuts)));
284  // Post the data
285  PostData(4,copyfCuts);
286  delete copyfCuts;
287  return;
288 }
289 //_________________________________________________
291 
292  //==================================================================================
293  // USER EXECUTION FUNCTION - start
294  //==================================================================================
295  //
296  // This is the main function for the heavy flavour analysis.
297  //
298  //==================================================================================
299 
300  if (!fInputEvent) {
301  Error("UserExec","NO EVENT FOUND!");
302  return;
303  }
304 
305  if(fEvents%100==0){
306  std::cout << "\r" << "Analysing event number: " << fEvents << std::endl;
307  }
308 
309  fEvents++;
310 
311  // Show trigger mask
312  if(fShowMask)
313  {
314  std::bitset<32> maskEV(((AliAODInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
315  std::cout << "Event mask: " << maskEV << std::endl;
316  std::cout << "Trigger mask: " << std::bitset<32>(fCuts->GetTriggerMask()) << std::endl;
317  }
318 
319  //==================================================================================
320  // EVENT INITIALIZATION - start
321  //==================================================================================
322 
323 
324  AliAODEvent * aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
325  TClonesArray * D0TracksFromFriendFile = 0;
326  fCEvents->Fill(1);
327 
328  if(!aodEvent && AODEvent() && IsStandardAOD())
329  {
330  // In case there is an AOD handler writing a standard AOD, use the AOD
331  // event in memory rather than the input (ESD) event.
332  aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
333  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
334  // have to taken from the AOD event hold by the AliAODExtension
335  AliAODHandler* aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
336  if(aodHandler->GetExtensions())
337  {
338  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
339  AliAODEvent *aodFromExt = ext->GetAOD();
340  D0TracksFromFriendFile=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
341  }
342  }
343  else
344  {
345  D0TracksFromFriendFile=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
346  }
347 
348  // fix for temporary bug in ESDfilter
349  // the AODs with null vertex pointer didn't pass the PhysSel
350  if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
351  fCEvents->Fill(2);
352 
353  fCounter->StoreEvent(aodEvent,fCuts,fUseMCInfo);
354 
355  // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
356  TString trigclass=aodEvent->GetFiredTriggerClasses();
357  if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")||trigclass.Contains("C0SMH-B-NOPF-ALL")) fCEvents->Fill(5);
358 
359  if(!fCuts->IsEventSelected(aodEvent))
360  {
361  if(fShowRejection) std::cout << "Event rejected by code: " << fCuts->GetWhyRejection() << std::endl;
362  if(fCuts->GetWhyRejection()==6) // rejected for Z vertex
363  {
364  fCEvents->Fill(6);
365  return;
366  }
367  }
368 
369  Bool_t isEvSel=fCuts->IsEventSelected(aodEvent);
370  if(!isEvSel) return;
371  fCEvents->Fill(3);
372 
373  //get the magnetic field
374  Double_t bz = (Double_t)aodEvent->GetMagneticField();
375 
376  // AOD primary vertex
377  AliAODVertex *primaryVertex = (AliAODVertex*)aodEvent->GetPrimaryVertex();
378  if(!primaryVertex) return;
379  if(primaryVertex->GetNContributors()<1) return;
380  fCEvents->Fill(4);
381 
382  if(!D0TracksFromFriendFile)
383  {
384  AliInfo("Could not find array of HF vertices, skipping the event");
385  return;
386  }
387  else AliDebug(2, Form("Found %d vertices",D0TracksFromFriendFile->GetEntriesFast()));
388 
389 
390  //==================================================================================
391  // EVENT INITIALIZATION - end
392  //==================================================================================
393  // B0 MC SIGNAL IDENTIFICATION - start
394  //==================================================================================
395 
396  // We create an array that contains all the monte carlo particles in the event
397  TClonesArray *mcTrackArray = 0;
398  if(fUseMCInfo) mcTrackArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
399 
400  // We create an array to save the MC labels of true signal tracks
401  TMatrix * B0toDStarPiLabelMatrix = new TMatrix(0,7);
402 
403  // We fill the array with all B0->DStarPi tracks
404  if(fUseMCInfo) {
405  B0toDStarPiSignalTracksInMC(mcTrackArray,aodEvent,B0toDStarPiLabelMatrix,fOutputB0MC);
406  }
407  //==================================================================================
408  // B0 MC SIGNAL IDENTIFICATION - end
409  //==================================================================================
410  // PARTICLE SELECTION LOOP - start
411  //==================================================================================
412  //
413  // Here we select and reconstruct the particles for the B0->D*Pion decay.
414  //
415  //==================================================================================
416 
417 
418  DStarPionSelection(aodEvent,primaryVertex,bz,mcTrackArray,B0toDStarPiLabelMatrix);
419  B0PionSelection(aodEvent,primaryVertex,bz,mcTrackArray,B0toDStarPiLabelMatrix);
420 
421  D0Selection(aodEvent,primaryVertex,bz,mcTrackArray,B0toDStarPiLabelMatrix,D0TracksFromFriendFile);
422  DStarAndB0Selection(aodEvent,primaryVertex,bz,mcTrackArray,B0toDStarPiLabelMatrix,D0TracksFromFriendFile);
423 
424  // Clear arrays and memory management:
425  fD0Tracks->erase(fD0Tracks->begin(),fD0Tracks->end());
426  fB0PionTracks->erase(fB0PionTracks->begin(),fB0PionTracks->end());
427  fDStarPionTracks->erase(fDStarPionTracks->begin(),fDStarPionTracks->end());
428 
429  delete B0toDStarPiLabelMatrix; B0toDStarPiLabelMatrix = NULL;
430 
431  //==================================================================================
432  // PARTICLE SELECTION LOOP - end
433  //==================================================================================
434 
435  PostData(1,fOutput);
436  PostData(5,fCounter);
437  PostData(6,fOutputD0Pion);
438  PostData(7,fOutputD0Kaon);
439  PostData(8,fOutputDStarPion);
440  PostData(9,fOutputB0Pion);
441  PostData(10,fOutputD0);
442  PostData(11,fOutputDStar);
443  PostData(12,fOutputB0);
444  PostData(13,fOutputD0_D0Pt);
445  PostData(14,fOutputD0_DStarPt);
446  PostData(15,fOutputDStar_DStarPt);
447  PostData(16,fOutputB0MC);
448 
449 
450  //==================================================================================
451  // USER EXECUTION FUNCTION - end
452  //==================================================================================
453 
454 }
455 //________________________________________ terminate ___________________________
460 
461  //Info("Terminate","");
462  AliAnalysisTaskSE::Terminate();
463 
464  fOutput = dynamic_cast<TList*> (GetOutputData(1));
465  if (!fOutput) {
466  printf("ERROR: fOutput not available\n");
467  return;
468  }
469 
470  fCEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fCEvents"));
471 
472  fOutputD0Pion = dynamic_cast<TList*> (GetOutputData(6));
473  if (!fOutputD0Pion) {
474  printf("ERROR: fOutputD0Pion not available\n");
475  return;
476  }
477  fOutputD0Kaon = dynamic_cast<TList*> (GetOutputData(7));
478  if (!fOutputD0Kaon) {
479  printf("ERROR: fOutputD0Kaon not available\n");
480  return;
481  }
482  fOutputDStarPion = dynamic_cast<TList*> (GetOutputData(8));
483  if (!fOutputDStarPion) {
484  printf("ERROR: fOutputDStarPion not available\n");
485  return;
486  }
487  fOutputB0Pion = dynamic_cast<TList*> (GetOutputData(9));
488  if (!fOutputB0Pion) {
489  printf("ERROR: fOutputB0Pion not available\n");
490  return;
491  }
492  fOutputD0 = dynamic_cast<TList*> (GetOutputData(10));
493  if (!fOutputD0) {
494  printf("ERROR: fOutputD0 not available\n");
495  return;
496  }
497  fOutputDStar = dynamic_cast<TList*> (GetOutputData(11));
498  if (!fOutputDStar) {
499  printf("ERROR: fOutputDStar not available\n");
500  return;
501  }
502  fOutputB0 = dynamic_cast<TList*> (GetOutputData(12));
503  if (!fOutputB0) {
504  printf("ERROR: fOutputB0 not available\n");
505  return;
506  }
507  fOutputD0_D0Pt = dynamic_cast<TList*> (GetOutputData(13));
508  if (!fOutputD0_D0Pt) {
509  printf("ERROR: fOutputD0_D0Pt not available\n");
510  return;
511  }
512  fOutputD0_DStarPt = dynamic_cast<TList*> (GetOutputData(14));
513  if (!fOutputD0_DStarPt) {
514  printf("ERROR: fOutputD0_DStarPt not available\n");
515  return;
516  }
517  fOutputDStar_DStarPt = dynamic_cast<TList*> (GetOutputData(15));
518  if (!fOutputDStar_DStarPt) {
519  printf("ERROR: fOutputDStar_DStarPt not available\n");
520  return;
521  }
522  fOutputB0MC = dynamic_cast<TList*> (GetOutputData(16));
523  if (!fOutputB0MC) {
524  printf("ERROR: fOutputB0MC not available\n");
525  return;
526  }
527  return;
528 }
529 
530 //___________________________________________________________________________
533  Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
534 
535  //slot #1
536  //OpenFile(1);
537  fOutput = new TList();
538  fOutput->SetOwner();
539  fOutput->SetName("chist0");
540 
541  fOutputD0Pion = new TList();
542  fOutputD0Pion->SetOwner();
543  fOutputD0Pion->SetName("listD0Pion");
544 
545  fOutputD0Kaon = new TList();
546  fOutputD0Kaon->SetOwner();
547  fOutputD0Kaon->SetName("listD0Kaon");
548 
549  fOutputDStarPion = new TList();
550  fOutputDStarPion->SetOwner();
551  fOutputDStarPion->SetName("listDStarPion");
552 
553  fOutputB0Pion = new TList();
554  fOutputB0Pion->SetOwner();
555  fOutputB0Pion->SetName("listB0Pion");
556 
557  fOutputD0 = new TList();
558  fOutputD0->SetOwner();
559  fOutputD0->SetName("listD0");
560 
561  fOutputDStar = new TList();
562  fOutputDStar->SetOwner();
563  fOutputDStar->SetName("listDStar");
564 
565  fOutputB0 = new TList();
566  fOutputB0->SetOwner();
567  fOutputB0->SetName("listB0");
568 
569  fOutputD0_D0Pt = new TList();
570  fOutputD0_D0Pt->SetOwner();
571  fOutputD0_D0Pt->SetName("listD0_D0Pt");
572 
573  fOutputD0_DStarPt = new TList();
574  fOutputD0_DStarPt->SetOwner();
575  fOutputD0_DStarPt->SetName("listD0_DStarPt");
576 
577  fOutputDStar_DStarPt = new TList();
578  fOutputDStar_DStarPt->SetOwner();
579  fOutputDStar_DStarPt->SetName("listDStar_DStarPt");
580 
581  fOutputB0MC = new TList();
582  fOutputB0MC->SetOwner();
583  fOutputB0MC->SetName("listB0MC");
584 
585  // define histograms
587 
588  //Counter for Normalization
589  fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));
590  fCounter->Init();
591 
592  PostData(1,fOutput);
593  PostData(6,fOutputD0Pion);
594  PostData(7,fOutputD0Kaon);
595  PostData(8,fOutputDStarPion);
596  PostData(9,fOutputB0Pion);
597  PostData(10,fOutputD0);
598  PostData(11,fOutputDStar);
599  PostData(12,fOutputB0);
600  PostData(13,fOutputD0_D0Pt);
601  PostData(14,fOutputD0_DStarPt);
602  PostData(15,fOutputDStar_DStarPt);
603  PostData(16,fOutputB0MC);
604 
605  return;
606 }
607 //___________________________________ histograms _______________________________________
610 
611 
612  fCEvents = new TH1F("fCEvents","conter",13,0,13);
613  fCEvents->SetStats(kTRUE);
614  fCEvents->GetXaxis()->SetTitle("1");
615  fCEvents->GetYaxis()->SetTitle("counts");
616  fCEvents->GetXaxis()->SetBinLabel(2,"no. of events");
617  fCEvents->GetXaxis()->SetBinLabel(3,"good prim vtx and B field");
618  fCEvents->GetXaxis()->SetBinLabel(4,"no event selected");
619  fCEvents->GetXaxis()->SetBinLabel(5,"no vtx contributors");
620  fCEvents->GetXaxis()->SetBinLabel(6,"trigger for PbPb");
621  fCEvents->GetXaxis()->SetBinLabel(7,"no z vtx");
622  fCEvents->GetXaxis()->SetBinLabel(12,"no. of D0 fail to be rec");
623  fOutput->Add(fCEvents);
624 
625  //====================================================
626 
627  TString name_mc_B0_pt ="mc_B0_pt";
628  TH1F* hist_mc_B0_pt = new TH1F(name_mc_B0_pt.Data(),"Pt monte carlo B0 in B0->D*#pi; p_{T} [GeV/c]; Entries",400,0,20);
629  hist_mc_B0_pt->Sumw2();
630  hist_mc_B0_pt->SetLineColor(6);
631  hist_mc_B0_pt->SetMarkerStyle(20);
632  hist_mc_B0_pt->SetMarkerSize(0.6);
633  hist_mc_B0_pt->SetMarkerColor(6);
634  TH1F* histogram_mc_B0_pt = (TH1F*)hist_mc_B0_pt->Clone();
635  fOutputB0MC->Add(histogram_mc_B0_pt);
636 
637  TString name_mc_B0_pion_pt ="mc_B0_pion_pt";
638  TH1F* hist_mc_B0_pion_pt = new TH1F(name_mc_B0_pion_pt.Data(),"Pt monte carlo pion of B0 in B0->D*#pi; p_{T} [GeV/c]; Entries",400,0,20);
639  hist_mc_B0_pion_pt->Sumw2();
640  hist_mc_B0_pion_pt->SetLineColor(6);
641  hist_mc_B0_pion_pt->SetMarkerStyle(20);
642  hist_mc_B0_pion_pt->SetMarkerSize(0.6);
643  hist_mc_B0_pion_pt->SetMarkerColor(6);
644  TH1F* histogram_mc_B0_pion_pt = (TH1F*)hist_mc_B0_pion_pt->Clone();
645  fOutputB0MC->Add(histogram_mc_B0_pion_pt);
646 
647  TString name_mc_DStar_pt ="mc_DStar_pt";
648  TH1F* hist_mc_DStar_pt = new TH1F(name_mc_DStar_pt.Data(),"Pt monte carlo DStar in B0->D*#pi; p_{T} [GeV/c]; Entries",400,0,20);
649  hist_mc_DStar_pt->Sumw2();
650  hist_mc_DStar_pt->SetLineColor(6);
651  hist_mc_DStar_pt->SetMarkerStyle(20);
652  hist_mc_DStar_pt->SetMarkerSize(0.6);
653  hist_mc_DStar_pt->SetMarkerColor(6);
654  TH1F* histogram_mc_DStar_pt = (TH1F*)hist_mc_DStar_pt->Clone();
655  fOutputB0MC->Add(histogram_mc_DStar_pt);
656 
657  TString name_mc_DStar_pion_pt ="mc_DStar_pion_pt";
658  TH1F* hist_mc_DStar_pion_pt = new TH1F(name_mc_DStar_pion_pt.Data(),"Pt monte carlo pion of DStar in B0->D*#pi; p_{T} [GeV/c]; Entries",400,0,20);
659  hist_mc_DStar_pion_pt->Sumw2();
660  hist_mc_DStar_pion_pt->SetLineColor(6);
661  hist_mc_DStar_pion_pt->SetMarkerStyle(20);
662  hist_mc_DStar_pion_pt->SetMarkerSize(0.6);
663  hist_mc_DStar_pion_pt->SetMarkerColor(6);
664  TH1F* histogram_mc_DStar_pion_pt = (TH1F*)hist_mc_DStar_pion_pt->Clone();
665  fOutputB0MC->Add(histogram_mc_DStar_pion_pt);
666 
667  TString name_mc_D0_pt ="mc_D0_pt";
668  TH1F* hist_mc_D0_pt = new TH1F(name_mc_D0_pt.Data(),"Pt monte carlo D0 in B0->D*#pi; p_{T} [GeV/c]; Entries",400,0,20);
669  hist_mc_D0_pt->Sumw2();
670  hist_mc_D0_pt->SetLineColor(6);
671  hist_mc_D0_pt->SetMarkerStyle(20);
672  hist_mc_D0_pt->SetMarkerSize(0.6);
673  hist_mc_D0_pt->SetMarkerColor(6);
674  TH1F* histogram_mc_D0_pt = (TH1F*)hist_mc_D0_pt->Clone();
675  fOutputB0MC->Add(histogram_mc_D0_pt);
676 
677  TString name_mc_D0_pion_pt ="mc_D0_pion_pt";
678  TH1F* hist_mc_D0_pion_pt = new TH1F(name_mc_D0_pion_pt.Data(),"Pt monte carlo pion of D0 in B0->D*#pi; p_{T} [GeV/c]; Entries",400,0,20);
679  hist_mc_D0_pion_pt->Sumw2();
680  hist_mc_D0_pion_pt->SetLineColor(6);
681  hist_mc_D0_pion_pt->SetMarkerStyle(20);
682  hist_mc_D0_pion_pt->SetMarkerSize(0.6);
683  hist_mc_D0_pion_pt->SetMarkerColor(6);
684  TH1F* histogram_mc_D0_pion_pt = (TH1F*)hist_mc_D0_pion_pt->Clone();
685  fOutputB0MC->Add(histogram_mc_D0_pion_pt);
686 
687  TString name_mc_D0_kaon_pt ="mc_D0_kaon_pt";
688  TH1F* hist_mc_D0_kaon_pt = new TH1F(name_mc_D0_kaon_pt.Data(),"Pt monte carlo kaon of D0 in B0->D*#pi; p_{T} [GeV/c]; Entries",400,0,20);
689  hist_mc_D0_kaon_pt->Sumw2();
690  hist_mc_D0_kaon_pt->SetLineColor(6);
691  hist_mc_D0_kaon_pt->SetMarkerStyle(20);
692  hist_mc_D0_kaon_pt->SetMarkerSize(0.6);
693  hist_mc_D0_kaon_pt->SetMarkerColor(6);
694  TH1F* histogram_mc_D0_kaon_pt = (TH1F*)hist_mc_D0_kaon_pt->Clone();
695  fOutputB0MC->Add(histogram_mc_D0_kaon_pt);
696 
697  //==================================================
698 
699  TString name_dca_D0_DStarPion ="dca_D0_DStarPion";
700  TH1F* hist_dca_D0_DStarPion = new TH1F(name_dca_D0_DStarPion.Data(),"dca_D0_DStarPion; DCA [cm]; Entries",1000,0,0.2);
701  hist_dca_D0_DStarPion->Sumw2();
702  hist_dca_D0_DStarPion->SetLineColor(6);
703  hist_dca_D0_DStarPion->SetMarkerStyle(20);
704  hist_dca_D0_DStarPion->SetMarkerSize(0.6);
705  hist_dca_D0_DStarPion->SetMarkerColor(6);
706  TH1F* histogram_dca_D0_DStarPion = (TH1F*)hist_dca_D0_DStarPion->Clone();
707  fOutputB0MC->Add(histogram_dca_D0_DStarPion);
708 
709  TString name_dca_D0_B0Pion ="dca_D0_B0Pion";
710  TH1F* hist_dca_D0_B0Pion = new TH1F(name_dca_D0_B0Pion.Data(),"dca_D0_B0Pion; DCA [cm]; Entries",1000,0,0.2);
711  hist_dca_D0_B0Pion->Sumw2();
712  hist_dca_D0_B0Pion->SetLineColor(6);
713  hist_dca_D0_B0Pion->SetMarkerStyle(20);
714  hist_dca_D0_B0Pion->SetMarkerSize(0.6);
715  hist_dca_D0_B0Pion->SetMarkerColor(6);
716  TH1F* histogram_dca_D0_B0Pion = (TH1F*)hist_dca_D0_B0Pion->Clone();
717  fOutputB0MC->Add(histogram_dca_D0_B0Pion);
718 
719  TString name_dca_DStarPion_B0Pion ="dca_DStarPion_B0Pion";
720  TH1F* hist_dca_DStarPion_B0Pion = new TH1F(name_dca_DStarPion_B0Pion.Data(),"dca_DStarPion_B0Pion; DCA [cm]; Entries",1000,0,0.2);
721  hist_dca_DStarPion_B0Pion->Sumw2();
722  hist_dca_DStarPion_B0Pion->SetLineColor(6);
723  hist_dca_DStarPion_B0Pion->SetMarkerStyle(20);
724  hist_dca_DStarPion_B0Pion->SetMarkerSize(0.6);
725  hist_dca_DStarPion_B0Pion->SetMarkerColor(6);
726  TH1F* histogram_dca_DStarPion_B0Pion = (TH1F*)hist_dca_DStarPion_B0Pion->Clone();
727  fOutputB0MC->Add(histogram_dca_DStarPion_B0Pion);
728 
729  TString name_dca_Signal_D0_DStarPion ="dca_Signal_D0_DStarPion";
730  TH1F* hist_dca_Signal_D0_DStarPion = new TH1F(name_dca_Signal_D0_DStarPion.Data(),"dca_Signal_D0_DStarPion; DCA [cm]; Entries",1000,0,0.2);
731  hist_dca_Signal_D0_DStarPion->Sumw2();
732  hist_dca_Signal_D0_DStarPion->SetLineColor(4);
733  hist_dca_Signal_D0_DStarPion->SetMarkerStyle(20);
734  hist_dca_Signal_D0_DStarPion->SetMarkerSize(0.6);
735  hist_dca_Signal_D0_DStarPion->SetMarkerColor(4);
736  TH1F* histogram_dca_Signal_D0_DStarPion = (TH1F*)hist_dca_Signal_D0_DStarPion->Clone();
737  fOutputB0MC->Add(histogram_dca_Signal_D0_DStarPion);
738 
739  TString name_dca_Signal_D0_B0Pion ="dca_Signal_D0_B0Pion";
740  TH1F* hist_dca_Signal_D0_B0Pion = new TH1F(name_dca_Signal_D0_B0Pion.Data(),"dca_Signal_D0_B0Pion; DCA [cm]; Entries",1000,0,0.2);
741  hist_dca_Signal_D0_B0Pion->Sumw2();
742  hist_dca_Signal_D0_B0Pion->SetLineColor(4);
743  hist_dca_Signal_D0_B0Pion->SetMarkerStyle(20);
744  hist_dca_Signal_D0_B0Pion->SetMarkerSize(0.6);
745  hist_dca_Signal_D0_B0Pion->SetMarkerColor(4);
746  TH1F* histogram_dca_Signal_D0_B0Pion = (TH1F*)hist_dca_Signal_D0_B0Pion->Clone();
747  fOutputB0MC->Add(histogram_dca_Signal_D0_B0Pion);
748 
749  TString name_dca_Signal_DStarPion_B0Pion ="dca_Signal_DStarPion_B0Pion";
750  TH1F* hist_dca_Signal_DStarPion_B0Pion = new TH1F(name_dca_Signal_DStarPion_B0Pion.Data(),"dca_Signal_DStarPion_B0Pion; DCA [cm]; Entries",1000,0,0.2);
751  hist_dca_Signal_DStarPion_B0Pion->Sumw2();
752  hist_dca_Signal_DStarPion_B0Pion->SetLineColor(4);
753  hist_dca_Signal_DStarPion_B0Pion->SetMarkerStyle(20);
754  hist_dca_Signal_DStarPion_B0Pion->SetMarkerSize(0.6);
755  hist_dca_Signal_DStarPion_B0Pion->SetMarkerColor(4);
756  TH1F* histogram_dca_Signal_DStarPion_B0Pion = (TH1F*)hist_dca_Signal_DStarPion_B0Pion->Clone();
757  fOutputB0MC->Add(histogram_dca_Signal_DStarPion_B0Pion);
758 
759  //==================================================
760 
761  TString name_B0s_in_analysis ="B0s_in_analysis";
762  TH1F* hist_B0s_in_analysis = new TH1F(name_B0s_in_analysis.Data(),"Number of B0 to kpipipi in the Analysis; Entries",10,0,10);
763  hist_B0s_in_analysis->Sumw2();
764  hist_B0s_in_analysis->SetLineColor(6);
765  hist_B0s_in_analysis->SetMarkerStyle(20);
766  hist_B0s_in_analysis->SetMarkerSize(0.6);
767  hist_B0s_in_analysis->SetMarkerColor(6);
768  hist_B0s_in_analysis->SetStats(kTRUE);
769  hist_B0s_in_analysis->GetXaxis()->SetBinLabel(1,"no. of B0s");
770  hist_B0s_in_analysis->GetXaxis()->SetBinLabel(2,"no. of B0s to kpipipi");
771  hist_B0s_in_analysis->GetXaxis()->SetBinLabel(3,"no. with all tracks in event");
772  hist_B0s_in_analysis->GetXaxis()->SetBinLabel(4,"no. ...");
773  hist_B0s_in_analysis->GetXaxis()->SetBinLabel(5,"no. ...");
774  hist_B0s_in_analysis->GetXaxis()->SetBinLabel(6,"no. ...");
775  hist_B0s_in_analysis->GetXaxis()->SetBinLabel(7,"no. ...");
776  hist_B0s_in_analysis->GetXaxis()->SetBinLabel(8,"no. ...");
777  TH1F* hist_B0s_in_analysis_mc = (TH1F*)hist_B0s_in_analysis->Clone();
778  fOutputB0MC->Add(hist_B0s_in_analysis_mc);
779 
780  TString name_B0s_per_bin ="B0s_per_bin";
781  TH1F* hist_B0s_per_bin = new TH1F(name_B0s_per_bin.Data(),"Number of B0 to kpipipi in the Analysis per bin; Entries",fnPtBins,0,fnPtBins);
782  for (Int_t i = 0; i < fnPtBins; ++i)
783  {
784  TString bin_name = "";
785  bin_name += fPtBinLimits[i];
786  bin_name += "-";
787  bin_name += fPtBinLimits[i+1];
788  hist_B0s_per_bin->GetXaxis()->SetBinLabel(i+1,bin_name);
789  }
790  TH1F* hist_B0s_per_bin_mc = (TH1F*)hist_B0s_per_bin->Clone();
791  fOutputB0MC->Add(hist_B0s_per_bin_mc);
792 
793  //======================================================================================================================================================
794 
795  //we make the histograms for the Pions and Kaon
796  for (Int_t i = 0; i < 4; i++){
797 
798  TString add_name = "";
799  TList * listout;
800  if(i==0) listout = fOutputD0Pion;
801  if(i==1) listout = fOutputD0Kaon;
802  if(i==2) listout = fOutputDStarPion;
803  if(i==3) listout = fOutputB0Pion;
804 
805  for (Int_t j = 0; j < 6; j++){
806  if(j==0) add_name = "";
807  if(j==1) add_name = "Signal";
808  if(j==2) add_name = "Cut";
809  if(j==3) add_name = "SignalCut";
810  if(j==4) add_name = "Result";
811  if(j==5) add_name = "SignalResult";
812 
813  TString name_Histogram = "";
814  TString discription_Histogram = "";
815  Int_t numberOfBins = 0;
816  Double_t lowerBound = 0.0;
817  Double_t upperBound = 0.0;
818 
819  for (Int_t k = 0; k < 9; ++k)
820  {
821  if(k==0){name_Histogram = "ptTrack"; discription_Histogram = "pt track; p_{T} [GeV/c]; Entries"; numberOfBins = 600; lowerBound = 0; upperBound = 30;}
822  if(k==1){name_Histogram = "momentumTrack"; discription_Histogram = "momentum track; p [GeV/c]; Entries"; numberOfBins = 600; lowerBound = 0; upperBound = 30;}
823  if(k==2){name_Histogram = "numberOfITS"; discription_Histogram = "Number of ITS clusters track; [#]; Entries"; numberOfBins = 10; lowerBound = -0.5; upperBound = 9.5;}
824  if(k==3){name_Histogram = "numberOfTPC"; discription_Histogram = "Number of TPC clusters track; [#]; Entries"; numberOfBins = 601; lowerBound = -0.5; upperBound = 600.5;}
825  if(k==4){name_Histogram = "pointsOnITS"; discription_Histogram = "Number of ITS clusters track per layer; [#]; Entries"; numberOfBins = 10; lowerBound = -0.5; upperBound = 9.5;}
826  if(k==5){name_Histogram = "nSigmaTPC"; discription_Histogram = "n sigma TPC for track PID; sigma; Entries"; numberOfBins = 500; lowerBound = -5; upperBound = 5;}
827  if(k==6){name_Histogram = "nSigmaTOF"; discription_Histogram = "n sigma TOF for track PID; sigma; Entries"; numberOfBins = 500; lowerBound = -5; upperBound = 5;}
828  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;}
829  if(k==8){name_Histogram = "impactParameter"; discription_Histogram = "Impact Parameter track; [cm]; Entries"; numberOfBins = 2000; lowerBound = 0; upperBound = 0.5;}
830 
831  name_Histogram += add_name;
832  TH1F* histogram = new TH1F(name_Histogram.Data(),discription_Histogram.Data(),numberOfBins,lowerBound,upperBound);
833  histogram->Sumw2();
834  if(j%2==0) histogram->SetLineColor(6);
835  if(j%2==1) histogram->SetLineColor(4);
836  histogram->SetMarkerStyle(20);
837  histogram->SetMarkerSize(0.6);
838  if(j%2==0) histogram->SetMarkerColor(6);
839  if(j%2==1) histogram->SetMarkerColor(4);
840  TH1F* histogram_Clone = (TH1F*)histogram->Clone();
841  listout->Add(histogram_Clone);
842  fDaughterHistogramArray[i][j][k] = (Int_t*)histogram_Clone;
843  }
844 
845  TString numberofparticlesperevent="numberofparticlesperevent";
846  numberofparticlesperevent += add_name;
847  TH1F* hist_numberofparticlesperevent = new TH1F(numberofparticlesperevent.Data(),"Number of particles per event; number of particles in one event; Entries",100,0,100);
848  hist_numberofparticlesperevent->Sumw2();
849  hist_numberofparticlesperevent->SetLineColor(6);
850  hist_numberofparticlesperevent->SetMarkerStyle(20);
851  hist_numberofparticlesperevent->SetMarkerSize(0.6);
852  hist_numberofparticlesperevent->SetMarkerColor(6);
853  TH1F* histogram_numberofparticlesperevent = (TH1F*)hist_numberofparticlesperevent->Clone();
854  listout->Add(histogram_numberofparticlesperevent);
855  fDaughterHistogramArray[i][j][12] = (Int_t*)histogram_numberofparticlesperevent;
856  }
857 
858  TH1F * effectOfCuts = new TH1F("effectOfCutsOnBackground","Removal counter",18,0,18);
859  effectOfCuts->SetStats(kTRUE);
860  effectOfCuts->GetXaxis()->SetTitle("Cut number");
861  effectOfCuts->GetYaxis()->SetTitle("Particles cut");
862  effectOfCuts->GetXaxis()->SetBinLabel(1,"total");
863  effectOfCuts->GetXaxis()->SetBinLabel(2,"1");
864  effectOfCuts->GetXaxis()->SetBinLabel(3,"2");
865  effectOfCuts->GetXaxis()->SetBinLabel(4,"3");
866  effectOfCuts->GetXaxis()->SetBinLabel(5,"4");
867  effectOfCuts->GetXaxis()->SetBinLabel(6,"5");
868  effectOfCuts->GetXaxis()->SetBinLabel(7,"6");
869  effectOfCuts->GetXaxis()->SetBinLabel(8,"7");
870  effectOfCuts->GetXaxis()->SetBinLabel(9,"8");
871  effectOfCuts->GetXaxis()->SetBinLabel(10,"9");
872  effectOfCuts->GetXaxis()->SetBinLabel(11,"10");
873  effectOfCuts->GetXaxis()->SetBinLabel(12,"11");
874  effectOfCuts->GetXaxis()->SetBinLabel(13,"12");
875  effectOfCuts->GetXaxis()->SetBinLabel(14,"13");
876  effectOfCuts->GetXaxis()->SetBinLabel(15,"14");
877  effectOfCuts->GetXaxis()->SetBinLabel(16,"15");
878  effectOfCuts->GetXaxis()->SetBinLabel(17,"16");
879  effectOfCuts->GetXaxis()->SetBinLabel(18,"17");
880  listout->Add(effectOfCuts);
881  fDaughterHistogramArray2D[i][0] = (Int_t*)effectOfCuts;
882 
883  TH1F * effectOfCutsMC = new TH1F("effectOfCutsOnSignal","Removal counter",18,0,18);
884  effectOfCutsMC->SetStats(kTRUE);
885  effectOfCutsMC->GetXaxis()->SetTitle("Cut number");
886  effectOfCutsMC->GetYaxis()->SetTitle("Particles cut");
887  effectOfCutsMC->GetXaxis()->SetBinLabel(1,"total");
888  effectOfCutsMC->GetXaxis()->SetBinLabel(2,"1");
889  effectOfCutsMC->GetXaxis()->SetBinLabel(3,"2");
890  effectOfCutsMC->GetXaxis()->SetBinLabel(4,"3");
891  effectOfCutsMC->GetXaxis()->SetBinLabel(5,"4");
892  effectOfCutsMC->GetXaxis()->SetBinLabel(6,"5");
893  effectOfCutsMC->GetXaxis()->SetBinLabel(7,"6");
894  effectOfCutsMC->GetXaxis()->SetBinLabel(8,"7");
895  effectOfCutsMC->GetXaxis()->SetBinLabel(9,"8");
896  effectOfCutsMC->GetXaxis()->SetBinLabel(10,"9");
897  effectOfCutsMC->GetXaxis()->SetBinLabel(11,"10");
898  effectOfCutsMC->GetXaxis()->SetBinLabel(12,"11");
899  effectOfCutsMC->GetXaxis()->SetBinLabel(13,"12");
900  effectOfCutsMC->GetXaxis()->SetBinLabel(14,"13");
901  effectOfCutsMC->GetXaxis()->SetBinLabel(15,"14");
902  effectOfCutsMC->GetXaxis()->SetBinLabel(16,"15");
903  effectOfCutsMC->GetXaxis()->SetBinLabel(17,"16");
904  effectOfCutsMC->GetXaxis()->SetBinLabel(18,"17");
905  listout->Add(effectOfCutsMC);
906  fDaughterHistogramArray2D[i][1] = (Int_t*)effectOfCutsMC;
907 
908  TString name_particle_pdg ="particle_pdg";
909  TH1F* hist_particle_pdg = new TH1F(name_particle_pdg.Data(),"Pdg code particle; pdg code; Entries",2000,-0.5,1999.5);
910  hist_particle_pdg->Sumw2();
911  hist_particle_pdg->SetLineColor(6);
912  hist_particle_pdg->SetMarkerStyle(20);
913  hist_particle_pdg->SetMarkerSize(0.6);
914  hist_particle_pdg->SetMarkerColor(6);
915  TH1F* histogram_particle_pdg = (TH1F*)hist_particle_pdg->Clone();
916  listout->Add(histogram_particle_pdg);
917  fDaughterHistogramArray2D[i][2] = (Int_t*)histogram_particle_pdg;
918 
919  TString name_particle_mother_pdg ="particle_mother_pdg";
920  TH1F* hist_particle_mother_pdg = new TH1F(name_particle_mother_pdg.Data(),"Pdg code particle mother; pdg code; Entries",2000,-0.5,1999.5);
921  hist_particle_mother_pdg->Sumw2();
922  hist_particle_mother_pdg->SetLineColor(6);
923  hist_particle_mother_pdg->SetMarkerStyle(20);
924  hist_particle_mother_pdg->SetMarkerSize(0.6);
925  hist_particle_mother_pdg->SetMarkerColor(6);
926  TH1F* histogram_particle_mother_pdg = (TH1F*)hist_particle_mother_pdg->Clone();
927  listout->Add(histogram_particle_mother_pdg);
928  fDaughterHistogramArray2D[i][3] = (Int_t*)histogram_particle_mother_pdg;
929 
930  TString name_ptB0_vs_ptTrack ="ptB0_vs_ptTrackBackground";
931  TH2F* hist_ptB0_vs_ptTrack = new TH2F(name_ptB0_vs_ptTrack.Data(),"Pt B0 vs Pt ; p_{T} B0 [GeV/c]; p_{T} track [GeV/c]",100,0,30,100,0,30);
932  hist_ptB0_vs_ptTrack->Sumw2();
933  hist_ptB0_vs_ptTrack->SetLineColor(6);
934  hist_ptB0_vs_ptTrack->SetMarkerStyle(20);
935  hist_ptB0_vs_ptTrack->SetMarkerSize(0.6);
936  hist_ptB0_vs_ptTrack->SetMarkerColor(6);
937  TH2F* histogram_ptB0_vs_ptTrack = (TH2F*)hist_ptB0_vs_ptTrack->Clone();
938  listout->Add(histogram_ptB0_vs_ptTrack);
939  fDaughterHistogramArray2D[i][4] = (Int_t*)histogram_ptB0_vs_ptTrack;
940 
941  TString name_ptB0_vs_ptTrackMC ="ptB0_vs_ptTrackSignal";
942  TH2F* hist_ptB0_vs_ptTrackMC = new TH2F(name_ptB0_vs_ptTrackMC.Data(),"Pt B0 vs Pt ; p_{T} B0 [GeV/c]; p_{T} track [GeV/c]",100,0,30,100,0,30);
943  hist_ptB0_vs_ptTrackMC->Sumw2();
944  hist_ptB0_vs_ptTrackMC->SetLineColor(4);
945  hist_ptB0_vs_ptTrackMC->SetMarkerStyle(20);
946  hist_ptB0_vs_ptTrackMC->SetMarkerSize(0.6);
947  hist_ptB0_vs_ptTrackMC->SetMarkerColor(6);
948  TH2F* histogram_ptB0_vs_ptTrackMC = (TH2F*)hist_ptB0_vs_ptTrackMC->Clone();
949  listout->Add(histogram_ptB0_vs_ptTrackMC);
950  fDaughterHistogramArray2D[i][5] = (Int_t*)histogram_ptB0_vs_ptTrackMC;
951  }
952 
953  //we make the histograms for the reconstructed particles
954  for (Int_t i = 0; i < 6; i++){
955 
956  TString add_name = "";
957  TList * listout;
958  Int_t nHistogramSets = 0;
959  if(i==0) {listout = fOutputD0; nHistogramSets = 6 + 2*fnPtBins;}
960  if(i==1) {listout = fOutputDStar; nHistogramSets = 6 + 2*fnPtBins;}
961  if(i==2) {listout = fOutputB0; nHistogramSets = 6 + 2*fnPtBins;}
962  if(i==3) {listout = fOutputD0_D0Pt; nHistogramSets = 2*fnPtBinsD0forD0ptbin;}
963  if(i==4) {listout = fOutputD0_DStarPt; nHistogramSets = 2*fnPtBinsD0forDStarptbin;}
964  if(i==5) {listout = fOutputDStar_DStarPt; nHistogramSets = 2*fnPtBinsDStarforDStarptbin;}
965 
966  for (Int_t j = 0; j < nHistogramSets; j++){
967  if(i<3)
968  {
969  if(j==0) add_name = "";
970  if(j==1) add_name = "Signal";
971  if(j==2) add_name = "Cut";
972  if(j==3) add_name = "SignalCut";
973  if(j==4) add_name = "Result";
974  if(j==5) add_name = "SignalResult";
975  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];}
976  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];}
977  }
978  if(i==3)
979  {
980  if(j%2==0) {add_name = "_ptbin_"; add_name += fPtBinLimitsD0forD0ptbin[j/2]; add_name += "_to_"; add_name += fPtBinLimitsD0forD0ptbin[1+j/2];}
981  if(j%2==1) {add_name = "Signal_ptbin_"; add_name += fPtBinLimitsD0forD0ptbin[(j-1)/2]; add_name += "_to_"; add_name += fPtBinLimitsD0forD0ptbin[1+(j-1)/2];}
982  }
983  if(i==4)
984  {
985  if(j%2==0) {add_name = "_ptbin_"; add_name += fPtBinLimitsD0forDStarptbin[j/2]; add_name += "_to_"; add_name += fPtBinLimitsD0forDStarptbin[1+j/2];}
986  if(j%2==1) {add_name = "Signal_ptbin_"; add_name += fPtBinLimitsD0forDStarptbin[(j-1)/2]; add_name += "_to_"; add_name += fPtBinLimitsD0forDStarptbin[1+(j-1)/2];}
987  }
988  if(i==5)
989  {
990  if(j%2==0) {add_name = "_ptbin_"; add_name += fPtBinLimitsDStarforDStarptbin[j/2]; add_name += "_to_"; add_name += fPtBinLimitsDStarforDStarptbin[1+j/2];}
991  if(j%2==1) {add_name = "Signal_ptbin_"; add_name += fPtBinLimitsDStarforDStarptbin[(j-1)/2]; add_name += "_to_"; add_name += fPtBinLimitsDStarforDStarptbin[1+(j-1)/2];}
992  }
993 
994 
995  TString name_Histogram = "";
996  TString discription_Histogram = "";
997  Int_t numberOfBins = 0;
998  Double_t lowerBound = 0.0;
999  Double_t upperBound = 0.0;
1000  Int_t numberOfBinsTwo = 0;
1001  Double_t lowerBoundTwo = 0.0;
1002  Double_t upperBoundTwo = 0.0;
1003 
1004  for (Int_t k = 0; k < 40; ++k)
1005  {
1006  if(k==0){name_Histogram = "ptMother"; discription_Histogram = "pt mother; p_{T} [GeV/c]; Entries"; numberOfBins = 300; lowerBound = 0; upperBound = 30;}
1007  if(k==1){name_Histogram = "ptFirstDaughter"; discription_Histogram = "pt first daughter; p_{T} [GeV/c]; Entries"; numberOfBins = 300; lowerBound = 0; upperBound = 30;}
1008  if(k==2){name_Histogram = "ptSecondDaughter"; discription_Histogram = "pt second daughter; p_{T} [GeV/c]; Entries"; numberOfBins = 300; lowerBound = 0; upperBound = 30;}
1009  if(k==3){name_Histogram = "etaMother"; discription_Histogram = "eta mother; #eta; Entries"; numberOfBins = 100; lowerBound = -2; upperBound = 2;}
1010  if(k==4){name_Histogram = "phiMother"; discription_Histogram = "phi mother; #phi; Entries"; numberOfBins = 25; lowerBound = 0; upperBound = 2*TMath::Pi();}
1011  if(k==5){name_Histogram = "d0Mother"; discription_Histogram = "d0 mother; [cm]; Entries"; numberOfBins = 2000; lowerBound = 0; upperBound = 0.5;}
1012  if(k==6){name_Histogram = "d0FirstDaughter"; discription_Histogram = "d0 first daughter; [cm]; Entries"; numberOfBins = 2000; lowerBound = 0; upperBound = 0.5;}
1013 
1014  if(k==7){name_Histogram = "d0SecondDaughter"; discription_Histogram = "d0 second daughter; [cm]; Entries"; numberOfBins = 2000; lowerBound = 0; upperBound = 0.5;}
1015 
1016  if(k==8){name_Histogram = "pointingAngleMother"; discription_Histogram = "pointing angle; [Cos(#theta)]; Entries"; numberOfBins = 100; lowerBound = -1; upperBound = 1;}
1017  if(k==9){name_Histogram = "impactProduct"; discription_Histogram = "impact product; [cm^{2}]; Entries"; numberOfBins = 500; lowerBound = -0.01; upperBound = 0.01;}
1018  if(k==10){name_Histogram = "impactProductXY"; discription_Histogram = "impact product XY; [cm^{2}]; Entries"; numberOfBins = 400; lowerBound = 0; upperBound = 0.1;}
1019  if(k==11){name_Histogram = "invariantMassMother"; discription_Histogram = "mass mother candidate; m [GeV/c^{2}]; Entries"; numberOfBins = 10000; lowerBound = 0; upperBound = 10;}
1020  if(k==12){name_Histogram = "deltaMassMother"; discription_Histogram = "mass mother candidate; m [GeV/c^{2}]; Entries"; numberOfBins = 10000; lowerBound = 0; upperBound = 10;}
1021  if(k==13){name_Histogram = "dcaMother"; discription_Histogram = "dca mother; distance [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 0.25;}
1022  if(k==14){name_Histogram = "vertexDistance"; discription_Histogram = "vertex distance between mother and primary vertex; distance [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 1;}
1023  if(k==15){name_Histogram = "normDecayLength"; discription_Histogram = "Normalized decay length w.r.t primary vertex; [cm]; Entries"; numberOfBins = 100; lowerBound = 0; upperBound = 50;}
1024  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;}
1025  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;}
1026  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;}
1027  if(k==19){name_Histogram = "angleMotherFirstDaughter"; discription_Histogram = "flight angle mother and first daughter; [Cos(#phi)]; Entries"; numberOfBins = 100; lowerBound = -1; upperBound = 1;}
1028  if(k==20){name_Histogram = "angleMotherSecondDaughter"; discription_Histogram = "flight angle mother and second daughter; [Cos(#phi)]; Entries"; numberOfBins = 100; lowerBound = 0.5; upperBound = 1;}
1029  if(k==21){name_Histogram = "angleBetweenBothDaughters"; discription_Histogram = "angle between both daughters; [Cos(#phi)]; Entries"; numberOfBins = 100; lowerBound = -1; upperBound = 1;}
1030  if(k==22){name_Histogram = "cosThetaStar"; discription_Histogram = "cosThetaStar; [Cos(#theta*)]; Entries"; numberOfBins = 200; lowerBound = -2; upperBound = 2;}
1031  if(k==23){name_Histogram = "vertexX"; discription_Histogram = "Vertex position; [cm]; Entries"; numberOfBins = 500; lowerBound = -5; upperBound = 5;}
1032  if(k==24){name_Histogram = "vertexY"; discription_Histogram = "Vertex position; [cm]; Entries"; numberOfBins = 500; lowerBound = -5; upperBound = 5;}
1033  if(k==25){name_Histogram = "vertexZ"; discription_Histogram = "Vertex position; [cm]; Entries"; numberOfBins = 500; lowerBound = -20; upperBound = 20;}
1034 
1035 
1036  if(k==26){if(i==0 || i==3 || i==4){name_Histogram = "pointingAngleToDStar"; discription_Histogram = "Pointing angle w.r.t. DStar decay vertex; [Cos(#theta)]; Entries"; numberOfBins = 200; lowerBound = -1; upperBound = 1;}
1037  else continue;}
1038  if(k==27){if(i==0 || i==3 || i==4){name_Histogram = "d0MotherToDStar"; discription_Histogram = "d0 Mother w.r.t. DStar decay vertex; [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 1;}
1039  else continue;}
1040  if(k==28){if(i==0 || i==3 || i==4){name_Histogram = "d0FirstDaughterToDStar"; discription_Histogram = "d0 first daughter w.r.t. DStar decay vertex; [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 1;}
1041  else continue;}
1042  if(k==29){if(i==0 || i==3 || i==4){name_Histogram = "d0SecondDaughterToDStar"; discription_Histogram = "d0 second daughter w.r.t. DStar decay vertex; [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 1;}
1043  else continue;}
1044  if(k==30){if(i==0 || i==3 || i==4){name_Histogram = "impactProductToDStar"; discription_Histogram = "impact product w.r.t. DStar decay vertex; [cm]; Entries"; numberOfBins = 500; lowerBound = -0.02; upperBound = 0.21;}
1045  else continue;}
1046  if(k==31){if(i==0 || i==3 || i==4){name_Histogram = "impactProductXYToDStar"; discription_Histogram = "impact product XY w.r.t. DStar decay vertex; [cm^{2}]; Entries"; numberOfBins = 100; lowerBound = 0; upperBound = 0.5;}
1047  else continue;}
1048  if(k==32){if(i==0 || i==3 || i==4){name_Histogram = "normDecayLengthToDStar"; discription_Histogram = "Normalized decay length w.r.t. DStar decay vertex; [cm]; Entries"; numberOfBins = 100; lowerBound = 0; upperBound = 50;}
1049  else continue;}
1050  if(k==33){if(i==0 || i==3 || i==4){name_Histogram = "pseudoProperDecayTimeToDStar"; discription_Histogram = "Pseudo Proper Decay Time w.r.t DStar vertex; [a.u.]; Entries"; numberOfBins = 1000; lowerBound = -1; upperBound = 1;}
1051  else continue;}
1052  if(k==34){if(i==0 || i==3 || i==4){name_Histogram = "DecayTimeToDStar"; discription_Histogram = "Decay Time w.r.t DStar vertex; [a.u.]; Entries"; numberOfBins = 100; lowerBound = 0; upperBound = 0.00000001;}
1053  else continue;}
1054  if(k==35){if(i==0 || i==3 || i==4){name_Histogram = "normDecayTimeToDStar"; discription_Histogram = "Normalized Decay Time w.r.t DStar vertex; [a.u.]; Entries"; numberOfBins = 100; lowerBound = 0; upperBound = 0.00001;}
1055  else continue;}
1056 
1057  if(k==36){name_Histogram = "topomaticFirstDaughter"; discription_Histogram = "topomatic d0 first daughter; [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 20;}
1058  if(k==37){name_Histogram = "topomaticSecondDaughter"; discription_Histogram = "topomatic d0 second daughter; [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 20;}
1059  if(k==38){name_Histogram = "topomaticMax"; discription_Histogram = "Max topomatic; [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 20;}
1060  if(k==39){name_Histogram = "topomaticMin"; discription_Histogram = "Min topomatic; [cm]; Entries"; numberOfBins = 500; lowerBound = 0; upperBound = 20;}
1061 
1062  name_Histogram += add_name;
1063  TH1F* histogram = new TH1F(name_Histogram.Data(),discription_Histogram.Data(),numberOfBins,lowerBound,upperBound);
1064  histogram->Sumw2();
1065  if(j%2==0) histogram->SetLineColor(6);
1066  if(j%2==1) histogram->SetLineColor(4);
1067  histogram->SetMarkerStyle(20);
1068  histogram->SetMarkerSize(0.6);
1069  if(j%2==0) histogram->SetMarkerColor(6);
1070  if(j%2==1) histogram->SetMarkerColor(4);
1071  TH1F* histogram_Clone = (TH1F*)histogram->Clone();
1072  listout->Add(histogram_Clone);
1073  fMotherHistogramArray[i][j][k] = (Int_t*)histogram_Clone;
1074  }
1075 
1076 
1077  name_Histogram = "";
1078  discription_Histogram = "";
1079  numberOfBins = 0;
1080  lowerBound = 0.0;
1081  upperBound = 0.0;
1082  numberOfBinsTwo = 0;
1083  lowerBoundTwo = 0.0;
1084  upperBoundTwo = 0.0;
1085 
1086  //we make the 2D histograms for the reconstructed particles
1087  Int_t nFirst = 0;
1088  Int_t nSecond = 1;
1089  Int_t nVariables = 8;
1090  Int_t nHistograms = nVariables * (nVariables - 1) / 2;
1091 
1092  TList * list2D = new TList();
1093  list2D->SetOwner();
1094  TString name2D = "2D_Histograms";
1095  name2D += add_name;
1096  list2D->SetName(name2D);
1097  listout->Add(list2D);
1098 
1099  for (Int_t k = 0; k < nHistograms; ++k)
1100  {
1101  numberOfBins = 50; numberOfBinsTwo = 50;
1102  if(nFirst==0){name_Histogram = "d0FirstDaughter"; discription_Histogram = "d0 first daughter [cm];"; lowerBound = 0; upperBound = 1;}
1103  if(nFirst==1){name_Histogram = "d0SecondDaughter"; discription_Histogram = "d0 second daughter [cm];"; lowerBound = 0; upperBound = 1;}
1104  if(nFirst==2){name_Histogram = "d0Mother"; discription_Histogram = "d0 mother [cm];"; lowerBound = 0; upperBound = 1;}
1105  if(nFirst==3){name_Histogram = "pointingAngleMother"; discription_Histogram = "pointing angle [Cos(#theta)];"; lowerBound = -1; upperBound = 1;}
1106  if(nFirst==4){name_Histogram = "impactProduct"; discription_Histogram = "impact product [cm^{2}];"; lowerBound = -0.01; upperBound = 0.01;}
1107  if(nFirst==5){name_Histogram = "impactProductXY"; discription_Histogram = "impact product XY [cm^{2}];"; lowerBound = 0; upperBound = 0.5;}
1108  if(nFirst==6){name_Histogram = "vertexDistance"; discription_Histogram = "vertex distance between mother and primary vertex [cm];"; lowerBound = 0; upperBound = 1;}
1109  if(nFirst==7){name_Histogram = "normDecayLength"; discription_Histogram = "Normalized decay length w.r.t primary vertex [cm];"; lowerBound = 0; upperBound = 50;}
1110  if(nFirst==8){name_Histogram = "pseudoProperDecayTime"; discription_Histogram = "Pseudo Proper Decay Time w.r.t primary vertex [a.u.];"; lowerBound = -10; upperBound = 10;}
1111 
1112  if(nSecond==0){name_Histogram += "d0FirstDaughter"; discription_Histogram += "d0 first daughter [cm];"; lowerBoundTwo = 0; upperBoundTwo = 1;}
1113  if(nSecond==1){name_Histogram += "d0SecondDaughter"; discription_Histogram += "d0 second daughter [cm];"; lowerBoundTwo = 0; upperBoundTwo = 1;}
1114  if(nSecond==2){name_Histogram += "d0Mother"; discription_Histogram += "d0 mother [cm];"; lowerBoundTwo = 0; upperBoundTwo = 1;}
1115  if(nSecond==3){name_Histogram += "pointingAngleMother"; discription_Histogram += "pointing angle [Cos(#theta)];"; lowerBoundTwo = -1; upperBoundTwo = 1;}
1116  if(nSecond==4){name_Histogram += "impactProduct"; discription_Histogram += "impact product [cm^{2}];"; lowerBoundTwo = -0.01; upperBoundTwo = 0.01;}
1117  if(nSecond==5){name_Histogram += "impactProductXY"; discription_Histogram += "impact product XY [cm^{2}];"; lowerBoundTwo = 0; upperBoundTwo = 0.5;}
1118  if(nSecond==6){name_Histogram += "vertexDistance"; discription_Histogram += "vertex distance between mother and primary vertex [cm];"; lowerBoundTwo = 0; upperBoundTwo = 1;}
1119  if(nSecond==7){name_Histogram += "normDecayLength"; discription_Histogram += "Normalized decay length w.r.t primary vertex [cm];"; lowerBoundTwo = 0; upperBoundTwo = 50;}
1120  if(nSecond==8){name_Histogram += "pseudoProperDecayTime"; discription_Histogram += "Pseudo Proper Decay Time w.r.t primary vertex [a.u.];"; lowerBoundTwo = -10; upperBoundTwo = 10;}
1121 
1122  name_Histogram += add_name;
1123  TH2F* histogram = new TH2F(name_Histogram.Data(),discription_Histogram.Data(),numberOfBins,lowerBound,upperBound,numberOfBinsTwo,lowerBoundTwo,upperBoundTwo);
1124  histogram->Sumw2();
1125  if(j%2==0) histogram->SetLineColor(6);
1126  if(j%2==1) histogram->SetLineColor(4);
1127  histogram->SetMarkerStyle(20);
1128  histogram->SetMarkerSize(0.6);
1129  histogram->SetMarkerColor(6);
1130  TH2F* histogram_Clone = (TH2F*)histogram->Clone();
1131  list2D->Add(histogram_Clone);
1132  fMotherHistogramArray2D[i][j][k] = (Int_t*)histogram_Clone;
1133 
1134  nSecond++;
1135  if(nSecond>nVariables)
1136  {
1137  nFirst++;
1138  nSecond = nFirst + 1;
1139  }
1140  }
1141  }
1142 
1143  TH1F * effectOfCuts = new TH1F("effectOfCuts","Removal counter",100,0,100);
1144  effectOfCuts->SetStats(kTRUE);
1145  effectOfCuts->GetXaxis()->SetTitle("Cut number");
1146  effectOfCuts->GetYaxis()->SetTitle("Particles cut");
1147  effectOfCuts->GetXaxis()->SetBinLabel(1,"total");
1148  for (Int_t i = 1; i < 100; ++i)
1149  {
1150  TString integerText = "";
1151  integerText += i;
1152  effectOfCuts->GetXaxis()->SetBinLabel(i+1,integerText);
1153  }
1154  listout->Add(effectOfCuts);
1155  fMotherHistogramArrayExtra[i][0] = (Int_t*)effectOfCuts;
1156 
1157  TH1F * effectOfCutsMC = new TH1F("effectOfCutsMC","Removal counter",100,0,100);
1158  effectOfCutsMC->SetStats(kTRUE);
1159  effectOfCutsMC->GetXaxis()->SetTitle("Cut number");
1160  effectOfCutsMC->GetYaxis()->SetTitle("Particles cut");
1161  effectOfCutsMC->GetXaxis()->SetBinLabel(1,"total");
1162  for (Int_t i = 1; i < 100; ++i)
1163  {
1164  TString integerText = "";
1165  integerText += i;
1166  effectOfCutsMC->GetXaxis()->SetBinLabel(i+1,integerText);
1167  }
1168  listout->Add(effectOfCutsMC);
1169  fMotherHistogramArrayExtra[i][1] = (Int_t*)effectOfCutsMC;
1170 
1171  TString name_particle_pdg ="particle_pdg";
1172  TH1F* hist_particle_pdg = new TH1F(name_particle_pdg.Data(),"Pdg code particle; pdg code; Entries",2000,-0.5,1999.5);
1173  hist_particle_pdg->Sumw2();
1174  hist_particle_pdg->SetLineColor(6);
1175  hist_particle_pdg->SetMarkerStyle(20);
1176  hist_particle_pdg->SetMarkerSize(0.6);
1177  hist_particle_pdg->SetMarkerColor(6);
1178  TH1F* histogram_particle_pdg = (TH1F*)hist_particle_pdg->Clone();
1179  listout->Add(histogram_particle_pdg);
1180  fMotherHistogramArrayExtra[i][2] = (Int_t*)histogram_particle_pdg;
1181 
1182  TString name_particle_mother_pdg ="particle_mother_pdg";
1183  TH1F* hist_particle_mother_pdg = new TH1F(name_particle_mother_pdg.Data(),"Pdg code particle mother; pdg code; Entries",2000,-0.5,1999.5);
1184  hist_particle_mother_pdg->Sumw2();
1185  hist_particle_mother_pdg->SetLineColor(6);
1186  hist_particle_mother_pdg->SetMarkerStyle(20);
1187  hist_particle_mother_pdg->SetMarkerSize(0.6);
1188  hist_particle_mother_pdg->SetMarkerColor(6);
1189  TH1F* histogram_particle_mother_pdg = (TH1F*)hist_particle_mother_pdg->Clone();
1190  listout->Add(histogram_particle_mother_pdg);
1191  fMotherHistogramArrayExtra[i][3] = (Int_t*)histogram_particle_mother_pdg;
1192 
1193  TString name_distance_vertex_from_real ="distance_vertex_from_real";
1194  TH1F* hist_distance_vertex_from_real = new TH1F(name_distance_vertex_from_real.Data(),"Distance reconstructed vertex from real vertex; distance [cm]; Entries",100,0,1);
1195  hist_distance_vertex_from_real->Sumw2();
1196  hist_distance_vertex_from_real->SetLineColor(6);
1197  hist_distance_vertex_from_real->SetMarkerStyle(20);
1198  hist_distance_vertex_from_real->SetMarkerSize(0.6);
1199  hist_distance_vertex_from_real->SetMarkerColor(6);
1200  TH1F* histogram_distance_vertex_from_real = (TH1F*)hist_distance_vertex_from_real->Clone();
1201  listout->Add(histogram_distance_vertex_from_real);
1202  fMotherHistogramArrayExtra[i][4] = (Int_t*)histogram_distance_vertex_from_real;
1203 
1204  TString name_distance_vertex_from_real_new ="distance_vertex_from_real_new";
1205  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",100,0,1);
1206  hist_distance_vertex_from_real_new->Sumw2();
1207  hist_distance_vertex_from_real_new->SetLineColor(6);
1208  hist_distance_vertex_from_real_new->SetMarkerStyle(20);
1209  hist_distance_vertex_from_real_new->SetMarkerSize(0.6);
1210  hist_distance_vertex_from_real_new->SetMarkerColor(6);
1211  TH1F* histogram_distance_vertex_from_real_new = (TH1F*)hist_distance_vertex_from_real_new->Clone();
1212  listout->Add(histogram_distance_vertex_from_real_new);
1213  fMotherHistogramArrayExtra[i][5] = (Int_t*)histogram_distance_vertex_from_real_new;
1214 
1215  TString name_momentum_resolution ="momentum_resolution";
1216  TH1F* hist_momentum_resolution = new TH1F(name_momentum_resolution.Data(),"Momentum resolution; difference between real and reconstructed momentum [GeV/c]; Entries",1000,0,1);
1217  hist_momentum_resolution->Sumw2();
1218  hist_momentum_resolution->SetLineColor(6);
1219  hist_momentum_resolution->SetMarkerStyle(20);
1220  hist_momentum_resolution->SetMarkerSize(0.6);
1221  hist_momentum_resolution->SetMarkerColor(6);
1222  TH1F* histogram_momentum_resolution = (TH1F*)hist_momentum_resolution->Clone();
1223  listout->Add(histogram_momentum_resolution);
1224  fMotherHistogramArrayExtra[i][6] = (Int_t*)histogram_momentum_resolution;
1225  }
1226 
1227  //we make the histograms for the same sign method histograms and the pt bins
1228  for (Int_t k = 0; k < fnPtBins+3; ++k){
1229  TString ptBinMother = "";
1230  if(k==0) ptBinMother = "";
1231  if(k==1) ptBinMother = "_ptbin_6_to_inf";
1232  if(k==2) ptBinMother = "_ptbin_3_to_inf";
1233  if(k>2) {ptBinMother += "_ptbin_"; ptBinMother += fPtBinLimits[k-3]; ptBinMother += "_to_"; ptBinMother += fPtBinLimits[k-2];}
1234 
1235  TString name_deltaMassMother ="deltaInvMassB0";
1236  name_deltaMassMother += ptBinMother;
1237  TH1F* hist_deltaMassMother = new TH1F(name_deltaMassMother.Data(),"mass mother candidate; m [GeV/c^2]; Entries",2000,0,20);
1238  hist_deltaMassMother->Sumw2();
1239  hist_deltaMassMother->SetLineColor(6);
1240  hist_deltaMassMother->SetMarkerStyle(20);
1241  hist_deltaMassMother->SetMarkerSize(0.6);
1242  hist_deltaMassMother->SetMarkerColor(6);
1243  TH1F* histogram_deltaMassMother = (TH1F*)hist_deltaMassMother->Clone();
1244  fOutputB0MC->Add(histogram_deltaMassMother);
1245 
1246  for (Int_t i = 0; i < 3; ++i){
1247  TString signName = "";
1248  if(i==0) signName = "";
1249  if(i==1) signName = "_SameSign";
1250  if(i==2) signName = "_SignSum";
1251  TString name_invariantMassMother ="invariantMassB0";
1252  name_invariantMassMother += ptBinMother + signName;
1253  TH1F* hist_invariantMassMother = new TH1F(name_invariantMassMother.Data(),"mass mother candidate; m [GeV/c^2]; Entries",2000,0,20);
1254  hist_invariantMassMother->Sumw2();
1255  hist_invariantMassMother->SetLineColor(6);
1256  hist_invariantMassMother->SetMarkerStyle(20);
1257  hist_invariantMassMother->SetMarkerSize(0.6);
1258  hist_invariantMassMother->SetMarkerColor(6);
1259  TH1F* histogram_invariantMassMother = (TH1F*)hist_invariantMassMother->Clone();
1260  fOutputB0MC->Add(histogram_invariantMassMother);
1261  }
1262  }
1263 
1264  TString name_cutEffectBackground ="cutEffectBackground";
1265  TH2I* hist_cutEffectBackground = new TH2I(name_cutEffectBackground.Data(),"Effect of Cuts on background; cut number; cut number",99,0,99,99,0,99);
1266  for (int i = 0; i < 99; ++i)
1267  {
1268  TString integerText = "";
1269  integerText += i;
1270  hist_cutEffectBackground->GetXaxis()->SetBinLabel(i+1,integerText);
1271  hist_cutEffectBackground->GetYaxis()->SetBinLabel(i+1,integerText);
1272  }
1273  TH2I* histogram_cutEffectBackground = (TH2I*)hist_cutEffectBackground->Clone();
1274  fOutputB0MC->Add(histogram_cutEffectBackground);
1275 
1276  TString name_cutEffectSignal ="cutEffectSignal";
1277  TH2I* hist_cutEffectSignal = new TH2I(name_cutEffectSignal.Data(),"Effect of Cuts on Signal; cut number; cut number",99,0,99,99,0,99);
1278  for (Int_t i = 0; i < 99; ++i)
1279  {
1280  TString integerText = "";
1281  integerText += i;
1282  hist_cutEffectSignal->GetXaxis()->SetBinLabel(i+1,integerText);
1283  hist_cutEffectSignal->GetYaxis()->SetBinLabel(i+1,integerText);
1284  }
1285  TH2I* histogram_cutEffectSignal = (TH2I*)hist_cutEffectSignal->Clone();
1286  fOutputB0MC->Add(histogram_cutEffectSignal);
1287 
1288  TString name_cutEffectUniqueBackground ="cutEffectUniqueBackground";
1289  TH1I* hist_cutEffectUniqueBackground = new TH1I(name_cutEffectUniqueBackground.Data(),"Effect of Cuts on Signal; cut number; cut number",99,0,99);
1290  for (Int_t i = 0; i < 99; ++i)
1291  {
1292  TString integerText = "";
1293  integerText += i;
1294  hist_cutEffectUniqueBackground->GetXaxis()->SetBinLabel(i+1,integerText);
1295  }
1296  TH1I* histogram_cutEffectUniqueBackground = (TH1I*)hist_cutEffectUniqueBackground->Clone();
1297  fOutputB0MC->Add(histogram_cutEffectUniqueBackground);
1298 
1299  TString name_cutEffectUniqueSignal ="cutEffectUniqueSignal";
1300  TH1I* hist_cutEffectUniqueSignal = new TH1I(name_cutEffectUniqueSignal.Data(),"Effect of Cuts on Signal; cut number; cut number",99,0,99);
1301  for (Int_t i = 0; i < 99; ++i)
1302  {
1303  TString integerText = "";
1304  integerText += i;
1305  hist_cutEffectUniqueSignal->GetXaxis()->SetBinLabel(i+1,integerText);
1306  }
1307  TH1I* histogram_cutEffectUniqueSignal = (TH1I*)hist_cutEffectUniqueSignal->Clone();
1308  fOutputB0MC->Add(histogram_cutEffectUniqueSignal);
1309 
1310  TString name_totalITSBackground ="totalITSBackground";
1311  TH1F* hist_totalITSBackground = new TH1F(name_totalITSBackground.Data(),"Total nr. of ITS hits for the daughters; number [#]; Entries",30,0,30);
1312  hist_totalITSBackground->Sumw2();
1313  hist_totalITSBackground->SetLineColor(6);
1314  hist_totalITSBackground->SetMarkerStyle(20);
1315  hist_totalITSBackground->SetMarkerSize(0.6);
1316  hist_totalITSBackground->SetMarkerColor(6);
1317  TH1F* histogram_totalITSBackground = (TH1F*)hist_totalITSBackground->Clone();
1318  fOutputB0MC->Add(histogram_totalITSBackground);
1319 
1320  TString name_totalITSSignal ="totalITSSignal";
1321  TH1F* hist_totalITSSignal = new TH1F(name_totalITSSignal.Data(),"Total nr. of ITS hits for the daughters; number [#]; Entries",30,0,30);
1322  hist_totalITSSignal->Sumw2();
1323  hist_totalITSSignal->SetLineColor(6);
1324  hist_totalITSSignal->SetMarkerStyle(20);
1325  hist_totalITSSignal->SetMarkerSize(0.6);
1326  hist_totalITSSignal->SetMarkerColor(6);
1327  TH1F* histogram_totalITSSignal = (TH1F*)hist_totalITSSignal->Clone();
1328  fOutputB0MC->Add(histogram_totalITSSignal);
1329 
1330  TString name_totalTPCBackground ="totalTPCBackground";
1331  TH1F* hist_totalTPCBackground = new TH1F(name_totalTPCBackground.Data(),"Total nr. of TPC hits for the daughters; number [#]; Entries",1000,0,1000);
1332  hist_totalTPCBackground->Sumw2();
1333  hist_totalTPCBackground->SetLineColor(6);
1334  hist_totalTPCBackground->SetMarkerStyle(20);
1335  hist_totalTPCBackground->SetMarkerSize(0.6);
1336  hist_totalTPCBackground->SetMarkerColor(6);
1337  TH1F* histogram_totalTPCBackground = (TH1F*)hist_totalTPCBackground->Clone();
1338  fOutputB0MC->Add(histogram_totalTPCBackground);
1339 
1340  TString name_totalTPCSignal ="totalTPCSignal";
1341  TH1F* hist_totalTPCSignal = new TH1F(name_totalTPCSignal.Data(),"Total nr. of TPC hits for the daughters; number [#]; Entries",1000,0,1000);
1342  hist_totalTPCSignal->Sumw2();
1343  hist_totalTPCSignal->SetLineColor(6);
1344  hist_totalTPCSignal->SetMarkerStyle(20);
1345  hist_totalTPCSignal->SetMarkerSize(0.6);
1346  hist_totalTPCSignal->SetMarkerColor(6);
1347  TH1F* histogram_totalTPCSignal = (TH1F*)hist_totalTPCSignal->Clone();
1348  fOutputB0MC->Add(histogram_totalTPCSignal);
1349 
1350  TString name_totalSigmaPIDBackground ="totalSigmaPIDBackground";
1351  TH1F* hist_totalSigmaPIDBackground = new TH1F(name_totalSigmaPIDBackground.Data(),"Total sigma of TPC and TOF PID for the daughters; number [#]; Entries",1000,0,100);
1352  hist_totalSigmaPIDBackground->Sumw2();
1353  hist_totalSigmaPIDBackground->SetLineColor(6);
1354  hist_totalSigmaPIDBackground->SetMarkerStyle(20);
1355  hist_totalSigmaPIDBackground->SetMarkerSize(0.6);
1356  hist_totalSigmaPIDBackground->SetMarkerColor(6);
1357  TH1F* histogram_totalSigmaPIDBackground = (TH1F*)hist_totalSigmaPIDBackground->Clone();
1358  fOutputB0MC->Add(histogram_totalSigmaPIDBackground);
1359 
1360  TString name_totalSigmaPIDSignal ="totalSigmaPIDSignal";
1361  TH1F* hist_totalSigmaPIDSignal = new TH1F(name_totalSigmaPIDSignal.Data(),"Total sigma of TPC and TOF PID for the daughters; number [#]; Entries",1000,0,100);
1362  hist_totalSigmaPIDSignal->Sumw2();
1363  hist_totalSigmaPIDSignal->SetLineColor(6);
1364  hist_totalSigmaPIDSignal->SetMarkerStyle(20);
1365  hist_totalSigmaPIDSignal->SetMarkerSize(0.6);
1366  hist_totalSigmaPIDSignal->SetMarkerColor(6);
1367  TH1F* histogram_totalSigmaPIDSignal = (TH1F*)hist_totalSigmaPIDSignal->Clone();
1368  fOutputB0MC->Add(histogram_totalSigmaPIDSignal);
1369  return;
1370 }
1371 //-------------------------------------------------------------------------------------
1372 AliAODVertex* AliAnalysisTaskSEB0toDStarPi::RecalculateVertex(const AliVVertex *primary,TObjArray *tracks,Double_t bField, Double_t dispersion){
1373  //
1374  // Helper function to recalculate a vertex.
1375  //
1376 
1377  AliESDVertex *vertexESD = 0;
1378  AliAODVertex *vertexAOD = 0;
1379 
1380  AliVertexerTracks *vertexer = new AliVertexerTracks(bField);
1381 
1382  vertexer->SetVtxStart((AliESDVertex*)primary); //primary vertex
1383  vertexESD = (AliESDVertex*)vertexer->VertexForSelectedESDTracks(tracks);
1384 
1385  delete vertexer; vertexer=NULL;
1386 
1387  if(!vertexESD) return vertexAOD;
1388 
1389 
1390  if(vertexESD->GetNContributors()!=tracks->GetEntriesFast())
1391  {
1392  delete vertexESD; vertexESD=NULL;
1393  return vertexAOD;
1394  }
1395 
1396  // convert to AliAODVertex
1397  Double_t pos[3],cov[6],chi2perNDF;
1398  for(Int_t a=0;a<3;a++)pos[a]=0.;
1399  for(Int_t b=0;b<6;b++)cov[b]=0.;
1400  chi2perNDF=0;
1401 
1402  vertexESD->GetXYZ(pos); // position
1403  vertexESD->GetCovMatrix(cov); //covariance matrix
1404 
1405 
1406  Double_t vertRadius2=pos[0]*pos[0]+pos[1]*pos[1];
1407  if(vertRadius2>8.) //(2.82)^2 radius beam pipe
1408  {
1409  delete vertexESD; vertexESD=NULL;
1410  return vertexAOD;
1411  }
1412 
1413  chi2perNDF = vertexESD->GetChi2toNDF();
1414  dispersion = vertexESD->GetDispersion();
1415  delete vertexESD; vertexESD=NULL;
1416  Int_t nprongs = 2; //tracks->GetEntriesFast();
1417  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
1418 
1419  return vertexAOD;
1420 }
1421 //-------------------------------------------------------------------------------------
1422 void AliAnalysisTaskSEB0toDStarPi::B0toDStarPiSignalTracksInMC(TClonesArray * mcTrackArray,AliAODEvent* aodevent,TMatrix * B0toDStarPiLabelMatrix, TList *listout){
1423 
1424  TMatrix &particleMatrix = *B0toDStarPiLabelMatrix;
1425  for (Int_t i=0; i<mcTrackArray->GetEntriesFast(); i++){
1426 
1427  Int_t mcLabelPionB0 = 0;
1428  Int_t mcLabelPionDStar = 0;
1429  Int_t mcLabelPionD0 = 0;
1430  Int_t mcLabelKaon = 0;
1431  Int_t mcLabelD0 = 0;
1432  Int_t mcLabelDStar = 0;
1433  Int_t mcLabelB0 = 0;
1434 
1435  Double_t ptMC[7] = {0.0};
1436 
1437 
1438  Bool_t mcPionB0Present = kFALSE;
1439  Bool_t mcPionDStarPresent = kFALSE;
1440  Bool_t mcPionD0Present = kFALSE;
1441  Bool_t mcKaonPresent = kFALSE;
1442 
1443  // Below, we find all the MC labels for the true signal tracks
1444  AliAODMCParticle *mcTrackParticle = dynamic_cast< AliAODMCParticle*>(mcTrackArray->At(i));
1445  Int_t pdgCodeMC=TMath::Abs(mcTrackParticle->GetPdgCode());
1446 
1447  if (pdgCodeMC==511){ //if the track is a B0 we look at its daughters
1448 
1449  mcLabelB0 = i;
1450  Int_t nDaughterB0 = mcTrackParticle->GetNDaughters();
1451  ptMC[0] = mcTrackParticle->Pt();
1452 
1453  TString fillthis= "B0s_in_analysis";
1454  ((TH1F*)(listout->FindObject(fillthis)))->Fill(0);
1455 
1456  if(nDaughterB0==2){
1457  for(Int_t iDaughterB0=0; iDaughterB0<2; iDaughterB0++){
1458 
1459  AliAODMCParticle* daughterB0 = (AliAODMCParticle*)mcTrackArray->At(mcTrackParticle->GetDaughter(iDaughterB0));
1460  if(!daughterB0) break;
1461  Int_t pdgCodeDaughterB0=TMath::Abs(daughterB0->GetPdgCode());
1462 
1463  if (pdgCodeDaughterB0==211){ //if the track is a pion we save its monte carlo label
1464  mcLabelPionB0 = mcTrackParticle->GetDaughter(iDaughterB0);
1465  mcPionB0Present = kTRUE;
1466  ptMC[1] = daughterB0->Pt();
1467 
1468  } else if (pdgCodeDaughterB0==413){ //if the track is a DStar we look at its daughters
1469  mcLabelDStar = mcTrackParticle->GetDaughter(iDaughterB0);
1470  Int_t nDaughterDStar = daughterB0->GetNDaughters();
1471  ptMC[2] = daughterB0->Pt();
1472 
1473  if(nDaughterDStar==2){
1474  for(Int_t iDaughterDStar=0; iDaughterDStar<2; iDaughterDStar++){
1475 
1476  AliAODMCParticle* daughterDStar = (AliAODMCParticle*)mcTrackArray->At(daughterB0->GetDaughter(iDaughterDStar));
1477  if(!daughterDStar) break;
1478  Int_t pdgCodeDaughterDStar=TMath::Abs(daughterDStar->GetPdgCode());
1479 
1480  if (pdgCodeDaughterDStar==211){ //if the track is a pion we save its monte carlo label
1481  mcLabelPionDStar = daughterB0->GetDaughter(iDaughterDStar);
1482  mcPionDStarPresent = kTRUE;
1483  ptMC[3] = daughterDStar->Pt();
1484 
1485  } else if (pdgCodeDaughterDStar==421){ //if the track is a D0 we look at its daughters
1486  mcLabelD0 = daughterB0->GetDaughter(iDaughterDStar);
1487  Int_t nDaughterD0 = daughterDStar->GetNDaughters();
1488  ptMC[4] = daughterDStar->Pt();
1489 
1490  if(nDaughterD0==2){
1491  for(Int_t iDaughterD0=0; iDaughterD0<2; iDaughterD0++){
1492 
1493  AliAODMCParticle* daughterD0 = (AliAODMCParticle*)mcTrackArray->At(daughterDStar->GetDaughter(iDaughterD0));
1494  if(!daughterD0) break;
1495  Int_t pdgCodeDaughterD0=TMath::Abs(daughterD0->GetPdgCode());
1496 
1497  if (pdgCodeDaughterD0==211){ //if the track is a pion we save its monte carlo label
1498  mcLabelPionD0 = daughterDStar->GetDaughter(iDaughterD0);
1499  ptMC[5] = daughterD0->Pt();
1500  mcPionD0Present = kTRUE;
1501 
1502  } else if (pdgCodeDaughterD0==321){ //if the track is a kaon we save its monte carlo label
1503  mcLabelKaon = daughterDStar->GetDaughter(iDaughterD0);;
1504  mcKaonPresent = kTRUE;
1505  ptMC[6] = daughterD0->Pt();
1506 
1507  } else break;
1508  }
1509  }
1510  } else break;
1511  }
1512  }
1513  } else break;
1514  }
1515  }
1516  }
1517  // Next, we save the labels to our array
1518  if (mcPionB0Present && mcPionDStarPresent && mcPionD0Present && mcKaonPresent){
1519  Int_t rows = B0toDStarPiLabelMatrix->GetNrows();
1520 
1521  B0toDStarPiLabelMatrix->ResizeTo(rows+1,7);
1522  particleMatrix(rows,0) = mcLabelPionB0;
1523  particleMatrix(rows,1) = mcLabelPionDStar;
1524  particleMatrix(rows,2) = mcLabelPionD0;
1525  particleMatrix(rows,3) = mcLabelKaon;
1526  particleMatrix(rows,4) = mcLabelD0;
1527  particleMatrix(rows,5) = mcLabelDStar;
1528  particleMatrix(rows,6) = mcLabelB0;
1529 
1530  // We also save information on the amount of signal tracks that exist in the MC dataset
1531  TString fillthis= "B0s_in_analysis";
1532  ((TH1F*)(listout->FindObject(fillthis)))->Fill(1);
1533 
1534  fillthis= "B0s_per_bin";
1535  for (Int_t j = 0; j < fnPtBins; ++j)
1536  {
1537  if(fPtBinLimits[j] < ptMC[0] && ptMC[0] < fPtBinLimits[j+1]) {((TH1F*)(listout->FindObject(fillthis)))->Fill(j); break;}
1538  }
1539 
1540  fillthis= "mc_B0_pt";
1541  ((TH1F*)(listout->FindObject(fillthis)))->Fill(ptMC[0]);
1542  fillthis= "mc_B0_pion_pt";
1543  ((TH1F*)(listout->FindObject(fillthis)))->Fill(ptMC[1]);
1544  fillthis= "mc_DStar_pt";
1545  ((TH1F*)(listout->FindObject(fillthis)))->Fill(ptMC[2]);
1546  fillthis= "mc_DStar_pion_pt";
1547  ((TH1F*)(listout->FindObject(fillthis)))->Fill(ptMC[3]);
1548  fillthis= "mc_D0_pt";
1549  ((TH1F*)(listout->FindObject(fillthis)))->Fill(ptMC[4]);
1550  fillthis= "mc_D0_pion_pt";
1551  ((TH1F*)(listout->FindObject(fillthis)))->Fill(ptMC[5]);
1552  fillthis= "mc_D0_kaon_pt";
1553  ((TH1F*)(listout->FindObject(fillthis)))->Fill(ptMC[6]);
1554 
1555  }
1556  }
1557 
1558  // Not all the tracks can be/are detected by the detector. We are only interested in tracks that are detectable and lie within the acceptance of the detector.
1559  // We remove the undetectable tracks from the array in order to get accurate information on the amount of signal that lies within the acceptance of our detector.
1560  Int_t numberOfB0s = 0;
1561  TArrayI correctLabelArray;
1562  for (Int_t i = 0; i < B0toDStarPiLabelMatrix->GetNrows(); i++)
1563  {
1564  Int_t particleCounter = 0;
1565  for (Int_t j = 0; j < 4; j++)
1566  {
1567  Int_t labelParticleInList = (Int_t)particleMatrix(i,j);
1568  for (Int_t k=0; k<aodevent->GetNumberOfTracks(); k++)
1569  {
1570  AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(aodevent->GetTrack(k));
1571  if(!aodTrack) AliFatal("Not a standard AOD");
1572  if(TMath::Abs(aodTrack->Eta())>0.8) continue;
1573  if(aodTrack->GetITSNcls() < 1) continue;
1574  if(aodTrack->GetTPCNcls() < 1) continue;
1575  if(aodTrack->GetLabel() == labelParticleInList)
1576  {
1577  particleCounter++;
1578  break;
1579  }
1580  }
1581  }
1582 
1583  if (particleCounter==4)
1584  {
1585  TString fillthis= "B0s_in_analysis";
1586  ((TH1F*)(listout->FindObject(fillthis)))->Fill(2);
1587  numberOfB0s++;
1588  correctLabelArray.Set(numberOfB0s);
1589  correctLabelArray.AddAt(i,numberOfB0s-1);
1590  }
1591  }
1592 
1593  for (Int_t i = 0; i < correctLabelArray.GetSize(); i++)
1594  {
1595  particleMatrix(i,0) = (Int_t)particleMatrix(correctLabelArray[i],0);
1596  particleMatrix(i,1) = (Int_t)particleMatrix(correctLabelArray[i],1);
1597  particleMatrix(i,2) = (Int_t)particleMatrix(correctLabelArray[i],2);
1598  particleMatrix(i,3) = (Int_t)particleMatrix(correctLabelArray[i],3);
1599  particleMatrix(i,4) = (Int_t)particleMatrix(correctLabelArray[i],4);
1600  particleMatrix(i,5) = (Int_t)particleMatrix(correctLabelArray[i],5);
1601  particleMatrix(i,6) = (Int_t)particleMatrix(correctLabelArray[i],6);
1602  }
1603  B0toDStarPiLabelMatrix->ResizeTo(correctLabelArray.GetSize(),7);
1604  return;
1605 }
1606 //-------------------------------------------------------------------------------------
1607 void AliAnalysisTaskSEB0toDStarPi::DStarPionSelection(AliAODEvent* aodEvent, AliAODVertex *primaryVertex, Double_t bz, TClonesArray * mcTrackArray, TMatrix * B0toDStarPiLabelMatrix){
1608 
1609  //we keep track of the number of particles we could use and how many we actually use after cuts
1610  Int_t numberofparticles = 0;
1611  Int_t numberofparticlesused = 0;
1612  Int_t iClonesArray = 0;
1613 
1614  TString fillthis = "";
1615 
1616  //we loop over all tracks in the event
1617  for (Int_t i=0; i<aodEvent->GetNumberOfTracks(); i++){
1618  AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(i));
1619  if(!aodTrack) AliFatal("Not a standard AOD");
1620 
1621  //quick quality cut
1622  if(aodTrack->GetITSNcls() < 1) continue;
1623  if(aodTrack->GetTPCNcls() < 1) continue;
1624  if(aodTrack->GetStatus()&AliESDtrack::kITSpureSA) continue;
1625  if(!(aodTrack->GetStatus()&AliESDtrack::kITSin)) continue;
1626  if(aodTrack->GetID() < 0) continue;
1627  Double_t covtest[21];
1628  if(!aodTrack->GetCovarianceXYZPxPyPz(covtest)) continue;
1629 
1630  Int_t mcLabelParticle = -1;
1631  Int_t pdgParticle = -1;
1632  mcLabelParticle = aodTrack->GetLabel();
1633 
1634  numberofparticles++;
1635 
1636  // we fill histograms with information of the track
1637  Double_t pt_track = aodTrack->Pt();
1638  Double_t momentum_track = aodTrack->P();
1639  Int_t numberOfITS = aodTrack->GetITSNcls();
1640  Int_t numberOfTPC = aodTrack->GetTPCNcls();
1641  Double_t nSigmaTPC = 0;
1642  Double_t nSigmaTOF = 0;
1643  Int_t pionPIDnumber = 2;
1644  Int_t kaonPIDnumber = 3;
1645  Int_t TPCok = 0;
1646  Int_t TOFok = 0;
1647 
1648  AliAODPidHF* trackPIDHF = (AliAODPidHF*)fCuts->GetPidHF();
1649  TPCok = trackPIDHF->GetnSigmaTPC(aodTrack, pionPIDnumber, nSigmaTPC);
1650  TOFok = trackPIDHF->GetnSigmaTOF(aodTrack, pionPIDnumber, nSigmaTOF);
1651 
1652  AliExternalTrackParam particleTrack;
1653  particleTrack.CopyFromVTrack(aodTrack);
1654  Double_t d0[2],covd0[3];
1655  particleTrack.PropagateToDCA(primaryVertex,bz,100.,d0,covd0);
1656 
1657  //we check if the particle is a signal track
1658  Bool_t isDesiredCandidate = kFALSE;
1659  if(fUseMCInfo){
1660  TMatrix &particleMatrix = *B0toDStarPiLabelMatrix;
1661  for (Int_t k = 0; k < B0toDStarPiLabelMatrix->GetNrows(); ++k){
1662  if(mcLabelParticle == (Int_t)particleMatrix(k,1)){
1663  isDesiredCandidate = kTRUE;
1664  break;
1665  }
1666  }
1667  }
1668 
1669  Int_t daughterType = 2;
1670 
1671 
1672  Int_t histType = 0;
1673  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
1674  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
1675  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
1676  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
1677 
1678  for (Int_t j = 0; j < 10; ++i)
1679  {
1680  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
1681 
1682  }
1683 
1684  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
1685  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
1686  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
1687  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
1688 
1689  if(isDesiredCandidate)
1690  {
1691  histType = 1;
1692  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
1693  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
1694  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
1695  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
1696 
1697  for (Int_t j = 0; j < 10; ++i)
1698  {
1699  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
1700 
1701  }
1702 
1703  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
1704  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
1705  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
1706  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
1707  }
1708 
1709  //we apply a number of cuts on the particle
1710  Bool_t bCut = kFALSE;
1711 
1712  //we apply a PID cut for a pion
1713  if(!(fCuts->SelectPID(aodTrack,pionPIDnumber))){
1714  if(isDesiredCandidate) {
1715  ((TH1F*)fDaughterHistogramArray2D[2][1])->Fill(2);
1716  } else ((TH1F*)fDaughterHistogramArray2D[2][0])->Fill(2);
1717  bCut = kTRUE;
1718  }
1719 
1720  if(aodTrack->GetITSNcls() < fCuts->GetMinITSNclsDStarPion()){
1721  if(isDesiredCandidate) {
1722  ((TH1F*)fDaughterHistogramArray2D[2][1])->Fill(3);
1723  } else ((TH1F*)fDaughterHistogramArray2D[2][0])->Fill(3);
1724  bCut = kTRUE;
1725  }
1726 
1727  if(aodTrack->GetTPCNcls() < fCuts->GetMinTPCNclsDStarPion()){
1728  if(isDesiredCandidate) {
1729  ((TH1F*)fDaughterHistogramArray2D[2][1])->Fill(4);
1730  } else ((TH1F*)fDaughterHistogramArray2D[2][0])->Fill(4);
1731  bCut = kTRUE;
1732  }
1733 
1734  if(fCuts->UseITSRefitDStarPion()==kTRUE){
1735  if(!(aodTrack->GetStatus()&AliESDtrack::kITSrefit)) {
1736  if(isDesiredCandidate) {
1737  ((TH1F*)fDaughterHistogramArray2D[2][1])->Fill(5);
1738  } else ((TH1F*)fDaughterHistogramArray2D[2][0])->Fill(5);
1739  bCut = kTRUE;
1740  }
1741  }
1742 
1743  if(fCuts->UseTPCRefitDStarPion()==kTRUE){
1744  if((!(aodTrack->GetStatus()&AliESDtrack::kTPCrefit))) {
1745  if(isDesiredCandidate) {
1746  ((TH1F*)fDaughterHistogramArray2D[2][1])->Fill(6);
1747  } else ((TH1F*)fDaughterHistogramArray2D[2][0])->Fill(6);
1748  bCut = kTRUE;
1749  }
1750  }
1751 
1752  if(fCuts->UseFilterBitDStarPion()==kTRUE){
1753  if(!(aodTrack->TestFilterMask(BIT(fCuts->GetFilterBitDStarPion())))) {
1754  if(isDesiredCandidate) {
1755  ((TH1F*)fDaughterHistogramArray2D[2][1])->Fill(7);
1756  } else ((TH1F*)fDaughterHistogramArray2D[2][0])->Fill(7);
1757  bCut = kTRUE;
1758  }
1759  }
1760 
1761  if(aodTrack->Pt() < fCuts->GetMinPtDStarPion()){
1762  if(isDesiredCandidate) {
1763  ((TH1F*)fDaughterHistogramArray2D[2][1])->Fill(8);
1764  } else ((TH1F*)fDaughterHistogramArray2D[2][0])->Fill(8);
1765  bCut = kTRUE;
1766  }
1767 
1768  if(!isDesiredCandidate && fQuickSignalAnalysis) bCut = kTRUE;
1769 
1770  if(bCut) {
1771  if(isDesiredCandidate) {
1772  ((TH1F*)fDaughterHistogramArray2D[2][1])->Fill(0);
1773  } else ((TH1F*)fDaughterHistogramArray2D[2][0])->Fill(0);
1774  continue;
1775  }
1776 
1777  //we fill histograms with track information of the tracks that pass the cuts
1778  histType = 2;
1779  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
1780  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
1781  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
1782  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
1783 
1784  for (Int_t j = 0; j < 10; ++i)
1785  {
1786  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
1787 
1788  }
1789 
1790  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
1791  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
1792  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
1793  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
1794 
1795  if(isDesiredCandidate)
1796  {
1797  histType = 3;
1798  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
1799  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
1800  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
1801  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
1802 
1803  for (Int_t j = 0; j < 10; ++i)
1804  {
1805  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
1806 
1807  }
1808 
1809  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
1810  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
1811  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
1812  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
1813  }
1814 
1815  fDStarPionTracks->push_back(i);
1816  numberofparticlesused++;
1817  }
1818 
1819  ((TH1F*)fDaughterHistogramArray[2][0][12])->Fill(numberofparticles);
1820  ((TH1F*)fDaughterHistogramArray[2][1][12])->Fill(numberofparticlesused);
1821  return;
1822 }
1823 //-------------------------------------------------------------------------------------
1824 void AliAnalysisTaskSEB0toDStarPi::B0PionSelection(AliAODEvent* aodEvent, AliAODVertex *primaryVertex, Double_t bz, TClonesArray * mcTrackArray, TMatrix * B0toDStarPiLabelMatrix){
1825 
1826  //we keep track of the number of particles we could use and how many we actually use after cuts
1827  Int_t numberofparticles = 0;
1828  Int_t numberofparticlesused = 0;
1829  Int_t iClonesArray = 0;
1830 
1831  TString fillthis = "";
1832 
1833  //we loop over all tracks in the event
1834  for (Int_t i=0; i<aodEvent->GetNumberOfTracks(); i++){
1835  AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(i));
1836  if(!aodTrack) AliFatal("Not a standard AOD");
1837 
1838  //quick quality cut
1839  if(aodTrack->GetITSNcls() < 1) continue;
1840  if(aodTrack->GetTPCNcls() < 1) continue;
1841  if(aodTrack->GetStatus()&AliESDtrack::kITSpureSA) continue;
1842  if(!(aodTrack->GetStatus()&AliESDtrack::kITSin)) continue;
1843  if(aodTrack->GetID() < 0) continue;
1844  Double_t covtest[21];
1845  if(!aodTrack->GetCovarianceXYZPxPyPz(covtest)) continue;
1846 
1847 
1848  Int_t mcLabelParticle = -1;
1849  Int_t pdgParticle = -1;
1850  mcLabelParticle = aodTrack->GetLabel();
1851 
1852  numberofparticles++;
1853 
1854  // we fill histograms with information of the track
1855  Double_t pt_track = aodTrack->Pt();
1856  Double_t momentum_track = aodTrack->P();
1857  Int_t numberOfITS = aodTrack->GetITSNcls();
1858  Int_t numberOfTPC = aodTrack->GetTPCNcls();
1859  Double_t nSigmaTPC = 0;
1860  Double_t nSigmaTOF = 0;
1861  Int_t pionPIDnumber = 2;
1862  Int_t kaonPIDnumber = 3;
1863  Int_t TPCok = 0;
1864  Int_t TOFok = 0;
1865 
1866  AliAODPidHF* trackPIDHF = (AliAODPidHF*)fCuts->GetPidHF();
1867  TPCok = trackPIDHF->GetnSigmaTPC(aodTrack, pionPIDnumber, nSigmaTPC);
1868  TOFok = trackPIDHF->GetnSigmaTOF(aodTrack, pionPIDnumber, nSigmaTOF);
1869 
1870  AliExternalTrackParam particleTrack;
1871  particleTrack.CopyFromVTrack(aodTrack);
1872  Double_t d0[2],covd0[3];
1873  particleTrack.PropagateToDCA(primaryVertex,bz,100.,d0,covd0);
1874 
1875 
1876  //we check if the particle is a signal track
1877  Bool_t isDesiredCandidate = kFALSE;
1878  if(fUseMCInfo){
1879  TMatrix &particleMatrix = *B0toDStarPiLabelMatrix;
1880  for (Int_t k = 0; k < B0toDStarPiLabelMatrix->GetNrows(); ++k){
1881  if(mcLabelParticle == (Int_t)particleMatrix(k,0)){
1882  isDesiredCandidate = kTRUE;
1883  break;
1884  }
1885  }
1886  }
1887 
1888  Int_t daughterType = 3;
1889 
1890  Int_t histType = 0;
1891  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
1892  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
1893  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
1894  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
1895 
1896  for (Int_t j = 0; j < 10; ++i)
1897  {
1898  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
1899 
1900  }
1901 
1902  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
1903  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
1904  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
1905  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
1906 
1907  if(isDesiredCandidate)
1908  {
1909  histType = 1;
1910  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
1911  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
1912  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
1913  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
1914 
1915  for (Int_t j = 0; j < 10; ++i)
1916  {
1917  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
1918 
1919  }
1920 
1921  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
1922  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
1923  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
1924  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
1925  }
1926 
1927  //we apply a number of cuts on the particle
1928  Bool_t bCut = kFALSE;
1929 
1930  //we apply a PID cut, 2 is used to indicate we look for a pion
1931  if(!(fCuts->SelectPID(aodTrack,2))){
1932  if(isDesiredCandidate) {
1933  ((TH1F*)fDaughterHistogramArray2D[3][1])->Fill(2);
1934  } else ((TH1F*)fDaughterHistogramArray2D[3][0])->Fill(2);
1935  bCut = kTRUE;
1936  }
1937 
1938  if(aodTrack->GetITSNcls() < fCuts->GetMinITSNclsB0Pion()){
1939  if(isDesiredCandidate) {
1940  ((TH1F*)fDaughterHistogramArray2D[3][1])->Fill(3);
1941  } else ((TH1F*)fDaughterHistogramArray2D[3][0])->Fill(3);
1942  bCut = kTRUE;
1943  }
1944 
1945  if(aodTrack->GetTPCNcls() < fCuts->GetMinTPCNclsB0Pion()){
1946  if(isDesiredCandidate) {
1947  ((TH1F*)fDaughterHistogramArray2D[3][1])->Fill(4);
1948  } else ((TH1F*)fDaughterHistogramArray2D[3][0])->Fill(4);
1949  bCut = kTRUE;
1950  }
1951 
1952  if(fCuts->UseITSRefitB0Pion()==kTRUE){
1953  if(!(aodTrack->GetStatus()&AliESDtrack::kITSrefit)) {
1954  if(isDesiredCandidate) {
1955  ((TH1F*)fDaughterHistogramArray2D[3][1])->Fill(5);
1956  } else ((TH1F*)fDaughterHistogramArray2D[3][0])->Fill(5);
1957  bCut = kTRUE;
1958  }
1959  }
1960 
1961  if(fCuts->UseTPCRefitB0Pion()==kTRUE){
1962  if((!(aodTrack->GetStatus()&AliESDtrack::kTPCrefit))) {
1963  if(isDesiredCandidate) {
1964  ((TH1F*)fDaughterHistogramArray2D[3][1])->Fill(6);
1965  } else ((TH1F*)fDaughterHistogramArray2D[3][0])->Fill(6);
1966  bCut = kTRUE;
1967  }
1968  }
1969 
1970  if(fCuts->UseFilterBitB0Pion()==kTRUE){
1971  if(!(aodTrack->TestFilterMask(BIT(fCuts->GetFilterBitB0Pion())))) {
1972  if(isDesiredCandidate) {
1973  ((TH1F*)fDaughterHistogramArray2D[3][1])->Fill(7);
1974  } else ((TH1F*)fDaughterHistogramArray2D[3][0])->Fill(7);
1975  bCut = kTRUE;
1976  }
1977  }
1978 
1979  if(aodTrack->Pt() < fCuts->GetMinPtB0Pion()){
1980  if(isDesiredCandidate) {
1981  ((TH1F*)fDaughterHistogramArray2D[3][1])->Fill(8);
1982  } else ((TH1F*)fDaughterHistogramArray2D[3][0])->Fill(8);
1983  bCut = kTRUE;
1984  }
1985 
1986  if(!isDesiredCandidate && fQuickSignalAnalysis) bCut = kTRUE;
1987 
1988  if(bCut) {
1989  if(isDesiredCandidate) {
1990  ((TH1F*)fDaughterHistogramArray2D[3][1])->Fill(0);
1991  } else ((TH1F*)fDaughterHistogramArray2D[3][0])->Fill(0);
1992  continue;
1993  }
1994 
1995  //we fill histograms with track information of the tracks that pass the cuts
1996  histType = 2;
1997  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
1998  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
1999  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
2000  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
2001 
2002  for (Int_t j = 0; j < 10; ++i)
2003  {
2004  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
2005 
2006  }
2007 
2008  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
2009  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
2010  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
2011  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
2012 
2013  if(isDesiredCandidate)
2014  {
2015  histType = 3;
2016  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
2017  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
2018  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
2019  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
2020 
2021  for (Int_t j = 0; j < 10; ++i)
2022  {
2023  if(aodTrack->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
2024 
2025  }
2026 
2027  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
2028  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
2029  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
2030  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0[0]);
2031  }
2032 
2033 
2034  fB0PionTracks->push_back(i);
2035  numberofparticlesused++;
2036  }
2037 
2038  ((TH1F*)fDaughterHistogramArray[3][0][12])->Fill(numberofparticles);
2039  ((TH1F*)fDaughterHistogramArray[3][1][12])->Fill(numberofparticlesused);
2040  return;
2041 }
2042 //-------------------------------------------------------------------------------------
2043 void AliAnalysisTaskSEB0toDStarPi::D0Selection(AliAODEvent* aodEvent, AliAODVertex *primaryVertex, Double_t bz,TClonesArray * mcTrackArray, TMatrix * B0toDStarPiLabelMatrix, TClonesArray * D0TracksFromFriendFile){
2044 
2045  TString fillthis = "";
2046 
2048 
2049  //next we loop over all the D0 candidates
2050  for (Int_t j = 0; j < D0TracksFromFriendFile->GetEntriesFast(); j++)
2051  {
2052 
2053  //we get the track of the D0
2054  AliAODRecoDecayHF2Prong * trackD0 = (AliAODRecoDecayHF2Prong*)(D0TracksFromFriendFile->At(j));
2055  if(!trackD0) {std::cout << "found none" << std::endl; continue;}
2056  if(trackD0 == NULL) {std::cout << "found NULL" << std::endl; continue;}
2057 
2058  if(!(vHF->FillRecoCand(aodEvent,trackD0))) //Fill the data members of the candidate only if they are empty.
2059  {
2060  fCEvents->Fill(12); //monitor how often this fails
2061  continue;
2062  }
2063 
2064  AliAODTrack * trackFirstDaughter = (AliAODTrack*)(trackD0->GetDaughter(0));
2065  AliAODTrack * trackSecondDaughter = (AliAODTrack*)(trackD0->GetDaughter(1));
2066 
2067  AliAODVertex *vertexMother = (AliAODVertex*)trackD0->GetSecondaryVtx();
2068 
2069  //we save the pdgcode of the used particle and its mother and check if it is a desired candidate
2070  Int_t pdgCodeMother = -1;
2071  Float_t pdgCodeGrandMother = -1;
2072  Bool_t isDesiredCandidate = kFALSE;
2073  Int_t motherType, histType;
2074  motherType = 0;
2075  Int_t mcLabelD0 = -1;
2076 
2077  if(fUseMCInfo)
2078  {
2079  mcLabelD0 = MatchCandidateToMonteCarlo(421,trackD0,mcTrackArray,B0toDStarPiLabelMatrix);
2080 
2081  if(mcLabelD0 >= 0)
2082  {
2083  isDesiredCandidate = kTRUE;
2084 
2085  Int_t mcLabelFirstTrack = -1;
2086  mcLabelFirstTrack = trackFirstDaughter->GetLabel();
2087 
2088  if(mcLabelFirstTrack >= 0)
2089  {
2090  AliAODMCParticle *mcParticleFirstTrack = (AliAODMCParticle*)mcTrackArray->At(mcLabelFirstTrack);
2091  AliAODMCParticle *mcMotherParticle = (AliAODMCParticle*)mcTrackArray->At(mcLabelD0);
2092 
2093  if(mcParticleFirstTrack && mcMotherParticle)
2094  {
2095  pdgCodeMother = mcMotherParticle->GetPdgCode();
2096 
2097  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()));
2098  ((TH1F*)fMotherHistogramArrayExtra[motherType][4])->Fill(vertex_distance);
2099 
2100  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()));
2101  ((TH1F*)fMotherHistogramArrayExtra[motherType][6])->Fill(momentum_resolution);
2102  }
2103  }
2104  }
2105  }
2106 
2107  // We fill the histograms
2108  histType = 0;
2109  FillD0Histograms(trackD0, primaryVertex, bz, motherType, histType);
2110  if(isDesiredCandidate && fUseMCInfo){
2111  histType = 1;
2112  FillD0Histograms(trackD0, primaryVertex, bz, motherType, histType,pdgCodeMother);
2113  }
2114 
2115  // Here we apply cuts on the particle
2116  Bool_t cutMother = kFALSE;
2117 
2118  Bool_t bCutArray[25] = {0};
2119  Int_t cutReturnValue = fCuts->IsD0forD0ptbinSelected(trackD0, 0, aodEvent, bCutArray);
2120  if(cutReturnValue == -1) cutMother = kTRUE;
2121  if(cutReturnValue == 0) cutMother = kTRUE;
2122 
2123 
2124  if(fGetCutInfo == kTRUE)
2125  {
2126  for (Int_t k = 0; k < 25; ++k)
2127  {
2128  if (bCutArray[k] == kTRUE){
2129  if(isDesiredCandidate){
2130  ((TH1F*)fMotherHistogramArray2D[motherType][1])->Fill(k+1);
2131  } else ((TH1F*)fMotherHistogramArray2D[motherType][0])->Fill(k+1);
2132  cutMother = kTRUE;
2133  }
2134  }
2135  }
2136 
2137 
2138  if(!isDesiredCandidate && fQuickSignalAnalysis) cutMother = kTRUE;
2139 
2140  if(cutMother){
2141  if(isDesiredCandidate){
2142  ((TH1F*)fMotherHistogramArray2D[motherType][1])->Fill(0);
2143  } else ((TH1F*)fMotherHistogramArray2D[motherType][0])->Fill(0);
2144  delete vertexMother; vertexMother = NULL;
2145  delete trackD0; trackD0 = NULL;
2146  continue;
2147  }
2148 
2149  // We fill the cut histograms
2150  histType = 2;
2151  FillD0Histograms(trackD0, primaryVertex, bz, motherType, histType);
2152  if(isDesiredCandidate && fUseMCInfo){
2153  histType = 3;
2154  FillD0Histograms(trackD0, primaryVertex, bz, motherType, histType,pdgCodeMother);
2155  }
2156 
2157  //we save the location of the D0 candidate
2158  fD0Tracks->push_back(j);
2159  }
2160 
2161  delete vHF;
2162  return;
2163 }
2164 //-------------------------------------------------------------------------------------
2165 void AliAnalysisTaskSEB0toDStarPi::DStarAndB0Selection(AliAODEvent* aodEvent, AliAODVertex *primaryVertex, Double_t bz,TClonesArray * mcTrackArray, TMatrix * B0toDStarPiLabelMatrix, TClonesArray * D0TracksFromFriendFile){
2166 
2167  TString fillthis = "";
2168 
2169  //we loop over all the DStar pion candidates
2170  for (Int_t i = 0; i < (Int_t)fDStarPionTracks->size(); i++)
2171  {
2172 
2173  //Save current Object count
2174  Int_t ObjectNumber = TProcessID::GetObjectCount();
2175 
2176  //we get the track of the DStar pion
2177  AliAODTrack * trackFirstDaughter = (AliAODTrack*)(aodEvent->GetTrack(fDStarPionTracks->at(i)));
2178  if(!trackFirstDaughter) continue;
2179 
2180  Int_t pdgD0 = 421;
2181  if(trackFirstDaughter->Charge() == -1) pdgD0 = -421;
2182 
2183  //next we loop over all the D0 candidates
2184  for (Int_t j = 0; j < (Int_t)fD0Tracks->size(); j++)
2185  {
2186 
2187  //we get the track of the D0
2188  AliAODRecoDecayHF2Prong * trackSecondDaughter = (AliAODRecoDecayHF2Prong*)(D0TracksFromFriendFile->At(fD0Tracks->at(j)));
2189  if(!trackSecondDaughter) {std::cout << "found none" << std::endl; continue;}
2190  if(trackSecondDaughter == NULL) {std::cout << "found NULL" << std::endl; continue;}
2191 
2192  //we check if the IDs of the tracks are different
2193  if(trackFirstDaughter->GetID() == trackSecondDaughter->GetProngID(0) || trackFirstDaughter->GetID() == trackSecondDaughter->GetProngID(1)) continue;
2194 
2195  //we check if the charges of the tracks are correct
2196  if(trackFirstDaughter->Charge() == trackSecondDaughter->Charge() || TMath::Abs(trackFirstDaughter->Charge() + trackSecondDaughter->Charge()) != 1) continue;
2197 
2198 
2199  //we check if the pions have the same charge
2200  if(trackFirstDaughter->Charge() == -1 && ((AliAODTrack*)trackSecondDaughter->GetDaughter(1))->Charge() != -1) continue;
2201  if(trackFirstDaughter->Charge() == 1 && ((AliAODTrack*)trackSecondDaughter->GetDaughter(0))->Charge() != 1) continue;
2202 
2203  //we apply a PID cut on the D0 daughters
2204 
2205  if(trackFirstDaughter->Charge()==1)
2206  {
2207  if(!(fCuts->SelectPID(((AliAODTrack*)trackSecondDaughter->GetDaughter(0)),2))) continue;
2208  if(!(fCuts->SelectPID(((AliAODTrack*)trackSecondDaughter->GetDaughter(1)),3))) continue;
2209  } else if (trackFirstDaughter->Charge()==-1){
2210  if(!(fCuts->SelectPID(((AliAODTrack*)trackSecondDaughter->GetDaughter(0)),3))) continue;
2211  if(!(fCuts->SelectPID(((AliAODTrack*)trackSecondDaughter->GetDaughter(1)),2))) continue;
2212  }
2213 
2214  //we make an estimate of the DStar vertex and make an initial broad invariant mass window cut
2215  AliExternalTrackParam DStarPionTrackParam;
2216  DStarPionTrackParam.CopyFromVTrack(trackFirstDaughter);
2217  AliExternalTrackParam D0TrackParam;
2218  D0TrackParam.CopyFromVTrack(trackSecondDaughter);
2219 
2220  // we calculate the vertex of the mother candidate
2221  TObjArray tracksTestVertex;
2222 
2223  tracksTestVertex.Add(&DStarPionTrackParam);
2224  tracksTestVertex.Add(&D0TrackParam);
2225 
2226  Double_t dispersionTest = 0;
2227  AliAODVertex *testVertex = RecalculateVertex(primaryVertex,&tracksTestVertex,bz,dispersionTest);
2228  if(!testVertex) {delete testVertex; testVertex = NULL; continue;}
2229 
2230  Double_t d0z0Test[2],covd0z0Test[3];
2231 
2232  //DStar creation with the new vertex
2233  DStarPionTrackParam.PropagateToDCA(testVertex,bz,100.,d0z0Test,covd0z0Test);
2234  D0TrackParam.PropagateToDCA(testVertex,bz,100.,d0z0Test,covd0z0Test);
2235  delete testVertex; testVertex = NULL;
2236 
2237  Double_t pdgMassPion = TDatabasePDG::Instance()->GetParticle(211)->Mass();
2238  Double_t pdgMassD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
2239  Double_t pdgMassDStar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
2240 
2241  Double_t energyDStarPion = pdgMassPion*pdgMassPion + DStarPionTrackParam.Px()*DStarPionTrackParam.Px()+DStarPionTrackParam.Py()*DStarPionTrackParam.Py()+DStarPionTrackParam.Pz()*DStarPionTrackParam.Pz();
2242  Double_t energyD0 = pdgMassD0*pdgMassD0 + D0TrackParam.Px()*D0TrackParam.Px()+D0TrackParam.Py()*D0TrackParam.Py()+D0TrackParam.Pz()*D0TrackParam.Pz();
2243  Double_t energySum = TMath::Sqrt(energyDStarPion) + TMath::Sqrt(energyD0);
2244 
2245  Double_t pxDStarTest = DStarPionTrackParam.Px() + D0TrackParam.Px();
2246  Double_t pyDStarTest = DStarPionTrackParam.Py() + D0TrackParam.Py();
2247  Double_t pzDStarTest = DStarPionTrackParam.Pz() + D0TrackParam.Pz();
2248  Double_t p2DStarTest = pxDStarTest*pxDStarTest + pyDStarTest*pyDStarTest + pzDStarTest*pzDStarTest;
2249 
2250  Double_t invMassDStarTest = TMath::Sqrt(energySum*energySum-p2DStarTest);
2251 
2252  //we use a mass window twice the size of the final cut.
2253  Int_t nCutIndex = 35;
2254  Bool_t bCutArrayTemp[25];
2255  Double_t cutVariableValue = TMath::Abs(invMassDStarTest-pdgMassDStar)/2.0;
2256  Bool_t bPassedCut = fCuts->ApplyCutOnVariableDStarforDStarptbin(nCutIndex,0,cutVariableValue,bCutArrayTemp);
2257  if(!bPassedCut) continue;
2258 
2259 
2260  //we loop over all the B0 pion candidates
2261  for (Int_t i = 0; i < (Int_t)fB0PionTracks->size(); i++)
2262  {
2263 
2264  //we get the track of the first daughter
2265  AliAODTrack * trackB0Pion = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(fB0PionTracks->at(i)));
2266  if(!trackB0Pion) continue;
2267 
2268  //we check if the IDs of the tracks are different
2269  AliAODTrack* twoProngdaughter0 = (AliAODTrack*)trackSecondDaughter->GetDaughter(0);
2270  AliAODTrack* twoProngdaughter1 = (AliAODTrack*)trackSecondDaughter->GetDaughter(1);
2271  UShort_t idProng0 = twoProngdaughter0->GetID();
2272  UShort_t idProng1 = twoProngdaughter1->GetID();
2273 
2274  if(trackB0Pion->GetID() == trackFirstDaughter->GetID() || trackB0Pion->GetID() == idProng0 || trackB0Pion->GetID() == idProng1) continue;
2275 
2276 
2277 
2278  //we check if the charges of the tracks are correct // later change this for like sign analysis.
2279  Bool_t bSameSign = kFALSE;
2280  if(trackB0Pion->Charge() == (trackSecondDaughter->Charge() + trackFirstDaughter->Charge()) && trackB0Pion->Charge() + (trackSecondDaughter->Charge() + trackFirstDaughter->Charge()) != 0) bSameSign = kTRUE;
2281 
2283  //
2284  // DStar Reconstruction
2285  //
2287 
2288  //we use the DStar pion, B0 pion, and D0 tracks to reconstruct the vertex for the B0 and DStar decay
2289  AliExternalTrackParam firstTrack;
2290  firstTrack.CopyFromVTrack(trackFirstDaughter);
2291  AliExternalTrackParam secondTrack;
2292  secondTrack.CopyFromVTrack(trackSecondDaughter);
2293  AliExternalTrackParam thirdTrack;
2294  thirdTrack.CopyFromVTrack(trackB0Pion);
2295 
2296  // we calculate the vertex
2297  TObjArray daughterTracksWithRecalculation;
2298 
2299  daughterTracksWithRecalculation.Add(&firstTrack);
2300  daughterTracksWithRecalculation.Add(&secondTrack);
2301  daughterTracksWithRecalculation.Add(&thirdTrack);
2302 
2303  Double_t dispersion = 0;
2304  AliAODVertex *vertexMother = RecalculateVertex(primaryVertex,&daughterTracksWithRecalculation,bz,dispersion);
2305  if(!vertexMother) {delete vertexMother; vertexMother = NULL; continue;}
2306 
2307 
2308  if(vertexMother->GetNDaughters()!=2)
2309  {
2310  std::cout << "bad reconstruction - number of daughters for vertex is incorrect" << std::endl;
2311  delete vertexMother; vertexMother = NULL;
2312  continue;
2313  }
2314 
2315  Double_t xdummyDStar=0.,ydummyDStar=0.,eDStar[2];
2316  Double_t d0z0DStar[2],covd0z0DStar[3],d0DStar[2],d0errDStar[2];
2317 
2318  //DStar creation with the new vertex
2319  firstTrack.PropagateToDCA(vertexMother,bz,100.,d0z0DStar,covd0z0DStar);
2320  secondTrack.PropagateToDCA(vertexMother,bz,100.,d0z0DStar,covd0z0DStar);
2321 
2322  Double_t pxDStar[2],pyDStar[2],pzDStar[2];
2323  pxDStar[0] = firstTrack.Px();
2324  pyDStar[0] = firstTrack.Py();
2325  pzDStar[0] = firstTrack.Pz();
2326  pxDStar[1] = secondTrack.Px();
2327  pyDStar[1] = secondTrack.Py();
2328  pzDStar[1] = secondTrack.Pz();
2329 
2330  Double_t xyz_track1[3];
2331  xyz_track1[0] = firstTrack.GetX();
2332  firstTrack.GetYAt(xyz_track1[0],bz,xyz_track1[1]);
2333  firstTrack.GetZAt(xyz_track1[0],bz,xyz_track1[2]);
2334 
2335  Double_t xyz_track2[3];
2336  xyz_track2[0] = secondTrack.GetX();
2337  secondTrack.GetYAt(xyz_track2[0],bz,xyz_track2[1]);
2338  secondTrack.GetZAt(xyz_track2[0],bz,xyz_track2[2]);
2339 
2340  Double_t distanceAtVertex = TMath::Sqrt((xyz_track1[0]-xyz_track2[0])*(xyz_track1[0]-xyz_track2[0]) + (xyz_track1[1]-xyz_track2[1])*(xyz_track1[1]-xyz_track2[1]) + (xyz_track1[2]-xyz_track2[2])*(xyz_track1[2]-xyz_track2[2]));
2341 
2342  firstTrack.PropagateToDCA(primaryVertex,bz,100.,d0z0DStar,covd0z0DStar);
2343  d0DStar[0] = d0z0DStar[0];
2344  d0errDStar[0] = TMath::Sqrt(covd0z0DStar[0]);
2345  secondTrack.PropagateToDCA(primaryVertex,bz,100.,d0z0DStar,covd0z0DStar);
2346  d0DStar[1] = d0z0DStar[0];
2347  d0errDStar[1] = TMath::Sqrt(covd0z0DStar[0]);
2348 
2349 
2350  Double_t dcaDStarPionD0 = secondTrack.GetDCA(&firstTrack,bz,xdummyDStar,ydummyDStar);
2351  Double_t dcaDStarPionB0Pion = secondTrack.GetDCA(&thirdTrack,bz,xdummyDStar,ydummyDStar);
2352  Double_t dcaB0PionD0 = thirdTrack.GetDCA(&firstTrack,bz,xdummyDStar,ydummyDStar);
2353 
2354  // if(dcaDStarPionD0 > 0.01 || dcaDStarPionB0Pion > 0.003 || dcaB0PionD0 > 0.01 )
2355  // {
2356  // delete vertexMother; vertexMother = NULL;
2357  // continue;
2358  // }
2359 
2360  Short_t chargeDStar = trackFirstDaughter->Charge() + trackSecondDaughter->Charge();
2361  AliAODVertex * vertexDStar = new AliAODVertex(*vertexMother);
2362  if(!vertexDStar)
2363  {
2364  std::cout << "no dstar vertex" << std::endl;
2365  delete vertexMother; vertexMother = NULL;
2366  delete vertexDStar; vertexDStar = NULL;
2367  continue;
2368 
2369  }
2370 
2371  Int_t nProngsDStar = 2;
2372  AliAODRecoDecayHF2Prong * trackDStar = new AliAODRecoDecayHF2Prong(vertexDStar,pxDStar,pyDStar,pzDStar,d0DStar,d0errDStar,distanceAtVertex);
2373  if(!trackDStar)
2374  {
2375  delete vertexMother; vertexMother = NULL;
2376  delete vertexDStar; vertexDStar = NULL;
2377  delete trackDStar; trackDStar = NULL;
2378  continue;
2379  }
2380 
2381  trackDStar->SetCharge(chargeDStar);
2382 
2383  UShort_t idDStar[2];
2384  idDStar[0]= trackFirstDaughter->GetID();
2385  idDStar[1]= 0;
2386 
2387  UInt_t prongsDStar[2];
2388  prongsDStar[0] = 211;
2389  prongsDStar[1] = 421;
2390 
2391  // cout << vertexDStar->GetNDaughters() << " ";
2392  if(vertexDStar->GetNDaughters()!=2)
2393  {
2394  std::cout << "bad reconstruction 2 - number of daughters for vertex is incorrect" << std::endl;
2395  delete vertexMother; vertexMother = NULL;
2396  delete vertexDStar; vertexDStar = NULL;
2397  delete trackDStar; trackDStar = NULL;
2398  continue;
2399  }
2400 
2401  // cout << trackDStar->GetSecondaryVtx() << " and " << ((AliAODVertex*)trackDStar->GetSecondaryVtx())->GetNDaughters() << " and " << trackFirstDaughter << endl;
2402 
2403  trackDStar->GetSecondaryVtx()->AddDaughter(trackFirstDaughter);
2404  trackDStar->GetSecondaryVtx()->AddDaughter(trackSecondDaughter);
2405  trackDStar->SetPrimaryVtxRef((AliAODVertex*)aodEvent->GetPrimaryVertex());
2406  trackDStar->SetProngIDs(2,idDStar);
2407 
2408 
2410  //
2411  // BO Reconstruction
2412  //
2414 
2415  // Use the new DStar candidate and the new vertex to create the B0 candidate
2416  Double_t xdummy=0.,ydummy=0.,dca,e[2];
2417  Double_t d0z0[2],covd0z0[3],d0[2],d0err[2];
2418 
2419  AliExternalTrackParam fourthTrack;
2420  fourthTrack.CopyFromVTrack(trackDStar);
2421 
2422  thirdTrack.PropagateToDCA(vertexMother,bz,100.,d0z0,covd0z0);
2423  fourthTrack.PropagateToDCA(vertexMother,bz,100.,d0z0,covd0z0);
2424 
2425  Double_t px[2],py[2],pz[2];
2426  px[0] = thirdTrack.Px();
2427  py[0] = thirdTrack.Py();
2428  pz[0] = thirdTrack.Pz();
2429  px[1] = fourthTrack.Px();
2430  py[1] = fourthTrack.Py();
2431  pz[1] = fourthTrack.Pz();
2432 
2433  UInt_t prongs[2];
2434  prongs[0] = 211;
2435  prongs[1] = 413;
2436 
2437  UShort_t id[2];
2438  id[0]= thirdTrack.GetID();
2439  id[1]= 0;
2440 
2441  thirdTrack.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
2442  d0[0] = d0z0[0];
2443  d0err[0] = TMath::Sqrt(covd0z0[0]);
2444  fourthTrack.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
2445  d0[1] = d0z0[0];
2446  d0err[1] = TMath::Sqrt(covd0z0[0]);
2447 
2448  dca = fourthTrack.GetDCA(&thirdTrack,bz,xdummy,ydummy);
2449 
2450 
2451  Short_t chargeMother = trackFirstDaughter->Charge() + trackDStar->Charge();
2452  Int_t nProngsB0 = 2;
2453  AliAODRecoDecayHF2Prong * trackB0 = new AliAODRecoDecayHF2Prong(vertexMother,px,py,pz,d0,d0err,dca);
2454  if(!trackB0)
2455  {
2456  delete vertexMother; vertexMother = NULL;
2457  delete vertexDStar; vertexDStar = NULL;
2458  delete trackDStar; trackDStar = NULL;
2459  delete trackB0; trackB0 = NULL;
2460  continue;
2461  }
2462 
2463  trackB0->SetCharge(chargeMother);
2464 
2465  trackB0->GetSecondaryVtx()->AddDaughter(trackB0Pion);
2466  trackB0->GetSecondaryVtx()->AddDaughter(trackDStar);
2467  trackB0->SetPrimaryVtxRef((AliAODVertex*)aodEvent->GetPrimaryVertex());
2468  trackB0->SetProngIDs(2,id);
2469 
2471  //
2472  // Cuts
2473  //
2475 
2476 
2477  // We check if the B0 candidate is a true signal in Monte Carlo
2478  Bool_t isDesiredCandidate = kFALSE;
2479  Int_t mcLabelB0 = -1;
2480  Int_t mcLabelDStar = -1;
2481  fillthis = "";
2482  Int_t motherType, histType;
2483  motherType = 1;
2484 
2485  if(fUseMCInfo)
2486  {
2487  mcLabelDStar = MatchCandidateToMonteCarlo(413,trackDStar,mcTrackArray,B0toDStarPiLabelMatrix);
2488  mcLabelB0 = MatchCandidateToMonteCarlo(511,trackB0,mcTrackArray,B0toDStarPiLabelMatrix);
2489 
2490  if (mcLabelB0 >= 0 && mcLabelDStar >= 0 && trackB0Pion->GetLabel() >= 0)
2491  {
2492  AliAODMCParticle *mcTrackDStarPion = (AliAODMCParticle*)mcTrackArray->At(trackB0Pion->GetLabel());
2493  AliAODMCParticle *mcTrackDStar = (AliAODMCParticle*)mcTrackArray->At(mcLabelDStar);
2494  // DStar
2495  Double_t vertex_distance = TMath::Sqrt((vertexMother->GetX() - mcTrackDStarPion->Xv())*(vertexMother->GetX() - mcTrackDStarPion->Xv()) + (vertexMother->GetY() - mcTrackDStarPion->Yv())*(vertexMother->GetY() - mcTrackDStarPion->Yv()) + (vertexMother->GetZ() - mcTrackDStarPion->Zv())*(vertexMother->GetZ() - mcTrackDStarPion->Zv()));
2496  ((TH1F*)fMotherHistogramArrayExtra[motherType][4])->Fill(vertex_distance);
2497 
2498  Double_t momentum_resolution = TMath::Sqrt((trackDStar->Px() - mcTrackDStar->Px())*(trackDStar->Px() - mcTrackDStar->Px()) + (trackDStar->Py() - mcTrackDStar->Py())*(trackDStar->Py() - mcTrackDStar->Py()) + (trackDStar->Pz() - mcTrackDStar->Pz())*(trackDStar->Pz() - mcTrackDStar->Pz()));
2499  ((TH1F*)fMotherHistogramArrayExtra[motherType][6])->Fill(momentum_resolution);
2500 
2501  isDesiredCandidate = kTRUE;
2502  }
2503  }
2504 
2505 
2506  // We fill the DStar histograms
2507  histType = 0;
2508  FillDStarAndB0Histograms(trackDStar, primaryVertex, bz, motherType, histType);
2509  if(isDesiredCandidate && fUseMCInfo)
2510  {
2511  histType = 1;
2512  FillDStarAndB0Histograms(trackDStar, primaryVertex, bz, motherType, histType);
2513  }
2514 
2515  // We apply cuts on the DStar
2516  Bool_t cutDStar = kFALSE;
2517 
2518  Bool_t bCutArrayDStar[25] = {0};
2519  Int_t cutReturnValueDStar = fCuts->IsDStarforDStarptbinSelected(trackDStar, 0, aodEvent, bCutArrayDStar);
2520  if(cutReturnValueDStar == -1) cutDStar = kTRUE;
2521  if(cutReturnValueDStar == 0) cutDStar = kTRUE;
2522 
2523  Bool_t bCutArrayD0[35] = {0};
2524  Int_t cutReturnValueD0 = fCuts->IsD0forDStarptbinSelected(trackDStar, 0, aodEvent, bCutArrayD0);
2525  if(cutReturnValueD0 == -1) cutDStar = kTRUE;
2526  if(cutReturnValueD0 == 0) cutDStar = kTRUE;
2527 
2528 
2529  if(fGetCutInfo == kTRUE)
2530  {
2531 
2532  for (Int_t n = 0; n < 25; ++n)
2533  {
2534  if(bCutArrayDStar[n] == kTRUE){
2535  if(isDesiredCandidate){
2536  ((TH1F*)fMotherHistogramArray2D[motherType][1])->Fill(n+1);
2537  } else ((TH1F*)fMotherHistogramArray2D[motherType][0])->Fill(n+1);
2538  cutDStar = kTRUE;
2539  }
2540  }
2541 
2542  for (Int_t n = 0; n < 35; ++n)
2543  {
2544  if(bCutArrayD0[n] == kTRUE){
2545  if(isDesiredCandidate){
2546  ((TH1F*)fMotherHistogramArray2D[motherType][1])->Fill(n+1+35);
2547  } else ((TH1F*)fMotherHistogramArray2D[motherType][0])->Fill(n+1+35);
2548  cutDStar = kTRUE;
2549  }
2550  }
2551  }
2552 
2553  if(!isDesiredCandidate && fQuickSignalAnalysis) cutDStar = kFALSE;
2554 
2555  if(cutDStar)
2556  {
2557  if(isDesiredCandidate)
2558  {
2559  ((TH1F*)fMotherHistogramArray2D[motherType][1])->Fill(0);
2560  } else ((TH1F*)fMotherHistogramArray2D[motherType][0])->Fill(0);
2561  delete vertexMother; vertexMother = NULL;
2562  delete vertexDStar; vertexDStar = NULL;
2563  delete trackDStar; trackDStar = NULL;
2564  continue;
2565  }
2566 
2567  // We fill the cut histograms
2568  histType = 2;
2569  FillDStarAndB0Histograms(trackDStar, primaryVertex, bz, motherType, histType);
2570  if(isDesiredCandidate && fUseMCInfo)
2571  {
2572  histType = 3;
2573  FillDStarAndB0Histograms(trackDStar, primaryVertex, bz, motherType, histType);
2574  }
2575 
2577  //
2578  // BO Reconstruction
2579  //
2580  //
2582 
2583 
2584  //we get information about the reconstructed B0
2585  Double_t ptMother = trackB0->Pt();
2586 
2587  fillthis = "";
2588 
2589  motherType = 2;
2590  histType = 0;
2591 
2592  if(isDesiredCandidate)
2593  {
2594  AliAODMCParticle *mcTrackFirstDaughter = (AliAODMCParticle*)mcTrackArray->At(trackB0Pion->GetLabel());
2595  AliAODMCParticle *mcTrackB0 = (AliAODMCParticle*)mcTrackArray->At(mcLabelB0);
2596 
2597  Double_t vertex_distance = TMath::Sqrt((vertexMother->GetX() - mcTrackFirstDaughter->Xv())*(vertexMother->GetX() - mcTrackFirstDaughter->Xv()) + (vertexMother->GetY() - mcTrackFirstDaughter->Yv())*(vertexMother->GetY() - mcTrackFirstDaughter->Yv()) + (vertexMother->GetZ() - mcTrackFirstDaughter->Zv())*(vertexMother->GetZ() - mcTrackFirstDaughter->Zv()));
2598  ((TH1F*)fMotherHistogramArrayExtra[motherType][4])->Fill(vertex_distance);
2599 
2600  Double_t momentum_resolution = TMath::Sqrt((trackB0->Px() - mcTrackB0->Px())*(trackB0->Px() - mcTrackB0->Px()) + (trackB0->Py() - mcTrackB0->Py())*(trackB0->Py() - mcTrackB0->Py()) + (trackB0->Pz() - mcTrackB0->Pz())*(trackB0->Pz() - mcTrackB0->Pz()));
2601  ((TH1F*)fMotherHistogramArrayExtra[motherType][6])->Fill(momentum_resolution);
2602  }
2603 
2604  if(!bSameSign)
2605  {
2606  // We fill the histograms
2607  histType = 0;
2608  FillDStarAndB0Histograms(trackB0, primaryVertex, bz, motherType, histType);
2609 
2610  if(isDesiredCandidate)
2611  {
2612  histType = 1;
2613  FillDStarAndB0Histograms(trackB0, primaryVertex, bz, motherType, histType);
2614  }
2615  }
2616 
2617 
2618  // We apply cuts
2619  Bool_t cutMother = kFALSE;
2620 
2621  Bool_t bCutArray[85] = {0};
2622  Int_t numberOfCuts = 85;
2623  Int_t cutReturnValue = fCuts->IsSelected(trackB0, 0, aodEvent, bCutArray);
2624  if(cutReturnValue == -1) cutMother = kTRUE;
2625  if(cutReturnValue == 0) cutMother = kTRUE;
2626 
2627 
2628  // We save information about the cuts
2629  TString histName = "";
2630  Double_t invariantMassMother = trackB0->InvMass(2,prongs);
2631  Double_t pdgMassMother=TDatabasePDG::Instance()->GetParticle(511)->Mass();
2632  Double_t massWindow = 0.125; //75; //GeV/c^2
2633  if(fGetCutInfo == kTRUE)
2634  {
2635  for (Int_t n = 0; n < 85; ++n)
2636  {
2637  if(bCutArray[n] == kTRUE){
2638  if(isDesiredCandidate){
2639  ((TH1F*)fMotherHistogramArray2D[motherType][1])->Fill(n+1);
2640  } else ((TH1F*)fMotherHistogramArray2D[motherType][0])->Fill(n+1);
2641  cutMother = kTRUE;
2642  }
2643  }
2644 
2645  if (TMath::Abs(invariantMassMother-pdgMassMother)<massWindow){
2646  for (Int_t i = 0; i < numberOfCuts; ++i) //total
2647  {
2648  if(bCutArray[i] == kFALSE) continue;
2649  for (Int_t j = 0; j < numberOfCuts; ++j)
2650  {
2651  if(bCutArray[j] == kFALSE) continue;
2652  if(isDesiredCandidate == kFALSE) histName ="cutEffectBackground";
2653  if(isDesiredCandidate == kTRUE) histName ="cutEffectSignal";
2654  ((TH2I*)(fOutputB0MC->FindObject(histName)))->Fill(i,j);
2655  }
2656  }
2657 
2658  for (Int_t i = 0; i < numberOfCuts; ++i) //unique
2659  {
2660  if(bCutArray[i] == kFALSE) continue;
2661  Bool_t bFill = kTRUE;
2662  for (Int_t j = 0; j < numberOfCuts; ++j)
2663  {
2664  if(i==j) continue;
2665  if(bCutArray[j] == kTRUE)
2666  {
2667  bFill = kFALSE;
2668  break;
2669  }
2670 
2671  }
2672  if(bFill == kTRUE)
2673  {
2674  if(isDesiredCandidate == kFALSE) histName ="cutEffectUniqueBackground";
2675  if(isDesiredCandidate == kTRUE) histName ="cutEffectUniqueSignal";
2676  ((TH1I*)(fOutputB0MC->FindObject(histName)))->Fill(i);
2677  }
2678  }
2679  }
2680  }
2681 
2682 
2683  if(!isDesiredCandidate && fQuickSignalAnalysis) cutMother = kTRUE;
2684 
2685  if(cutMother)
2686  {
2687  if(isDesiredCandidate)
2688  {
2689  ((TH1F*)fMotherHistogramArray2D[motherType][1])->Fill(0);
2690  } else ((TH1F*)fMotherHistogramArray2D[motherType][0])->Fill(0);
2691  delete vertexMother; vertexMother = NULL;
2692  delete trackB0; trackB0 = NULL;
2693  delete vertexDStar; vertexDStar = NULL;
2694  delete trackDStar; trackDStar = NULL;
2695  continue;
2696  }
2697 
2698  // We save the DCA information
2699  TString name_dca_D0_DStarPion ="dca_D0_DStarPion";
2700  ((TH1F*)(fOutputB0MC->FindObject(name_dca_D0_DStarPion)))->Fill(dcaDStarPionD0);
2701 
2702  TString name_dca_D0_B0Pion ="dca_D0_B0Pion";
2703  ((TH1F*)(fOutputB0MC->FindObject(name_dca_D0_B0Pion)))->Fill(dcaB0PionD0);
2704 
2705  TString name_dca_DStarPion_B0Pion ="dca_DStarPion_B0Pion";
2706  ((TH1F*)(fOutputB0MC->FindObject(name_dca_DStarPion_B0Pion)))->Fill(dcaDStarPionB0Pion);
2707 
2708  if(isDesiredCandidate && fUseMCInfo)
2709  {
2710  TString name_dca_MC_D0_DStarPion ="dca_MC_D0_DStarPion";
2711  ((TH1F*)(fOutputB0MC->FindObject(name_dca_MC_D0_DStarPion)))->Fill(dcaDStarPionD0);
2712 
2713  TString name_dca_MC_D0_B0Pion ="dca_MC_D0_B0Pion";
2714  ((TH1F*)(fOutputB0MC->FindObject(name_dca_MC_D0_B0Pion)))->Fill(dcaB0PionD0);
2715 
2716  TString name_dca_MC_DStarPion_B0Pion ="dca_MC_DStarPion_B0Pion";
2717  ((TH1F*)(fOutputB0MC->FindObject(name_dca_MC_DStarPion_B0Pion)))->Fill(dcaDStarPionB0Pion);
2718  }
2719 
2720 
2721  if(!bSameSign)
2722  {
2723  histType = 2;
2724  FillDStarAndB0Histograms(trackB0, primaryVertex, bz, motherType, histType);
2725  if(fUseMCInfo && isDesiredCandidate)
2726  {
2727  //fill mc histograms
2728  histType = 3;
2729  FillDStarAndB0Histograms(trackB0, primaryVertex, bz, motherType, histType);
2730  }
2731  }
2732 
2733  AliAODRecoDecayHF2Prong* selectedDStar = (AliAODRecoDecayHF2Prong*)trackB0->GetDaughter(1);
2734  AliAODRecoDecayHF2Prong* selectedD0 = (AliAODRecoDecayHF2Prong*)selectedDStar->GetDaughter(1);
2735 
2736  // We fill the final cut histograms with the candidates that have an invariant mass close to the PDG value. This way the background we use for optimizing the cuts will not be contaminated with candidates that don't affect the signal region.
2737  massWindow = 0.125; // GeV/c^2
2738  if (TMath::Abs(invariantMassMother-pdgMassMother)<massWindow)
2739  {
2740  if(!bSameSign)
2741  {
2742  FillFinalTrackHistograms(trackB0, isDesiredCandidate, mcTrackArray);
2743  if(!isDesiredCandidate)
2744  {
2745  motherType = 0; histType = 4; FillD0Histograms(selectedD0, primaryVertex, bz, motherType, histType, pdgD0);
2746  motherType = 1; histType = 4; FillDStarAndB0Histograms(selectedDStar, primaryVertex, bz, motherType, histType);
2747  motherType = 2; histType = 4; FillDStarAndB0Histograms(trackB0, primaryVertex, bz, motherType, histType);
2748  }
2749  if(isDesiredCandidate)
2750  {
2751  motherType = 0; histType = 5; FillD0Histograms(selectedD0, primaryVertex, bz, motherType, histType, pdgD0);
2752  motherType = 1; histType = 5; FillDStarAndB0Histograms(selectedDStar, primaryVertex, bz, motherType, histType);
2753  motherType = 2; histType = 5; FillDStarAndB0Histograms(trackB0, primaryVertex, bz, motherType, histType);
2754  }
2755  }
2756  }
2757 
2758  //Here we fill the histograms per pt bin and apply the same sign method
2759  TString ptBinMother = "";
2760  Int_t ptBin = fCuts->PtBin(trackB0->Pt());
2761  ptBinMother += "_ptbin_"; ptBinMother += fPtBinLimits[ptBin]; ptBinMother += "_to_"; ptBinMother += fPtBinLimits[ptBin+1];
2762  histType = 6 + 2 * ptBin;
2763 
2764  Int_t d0PtBin = fCuts->PtBinD0forD0ptbin(selectedD0->Pt());
2765  Int_t histTypeD0 = 2 * d0PtBin;
2766 
2767  Int_t d0DStarPtBin = fCuts->PtBinD0forDStarptbin(selectedDStar->Pt());
2768  Int_t histTypeD0DStar = 2 * d0DStarPtBin;
2769 
2770  Int_t dstarPtBin = fCuts->PtBinDStarforDStarptbin(selectedDStar->Pt());
2771  Int_t histTypeDStar = 2 * dstarPtBin;
2772 
2773 
2774  if (TMath::Abs(invariantMassMother-pdgMassMother)<massWindow)
2775  {
2776  if(!bSameSign && histType > 5)
2777  {
2778  if(!isDesiredCandidate)
2779  {
2780  motherType = 0; FillD0Histograms(selectedD0, primaryVertex, bz, motherType, histType, pdgD0);
2781  motherType = 1; FillDStarAndB0Histograms(selectedDStar, primaryVertex, bz, motherType, histType);
2782  motherType = 2; FillDStarAndB0Histograms(trackB0, primaryVertex, bz, motherType, histType);
2783  motherType = 3; FillD0Histograms(selectedD0, primaryVertex, bz, motherType, histTypeD0, pdgD0);
2784  motherType = 4; FillD0Histograms(selectedD0, primaryVertex, bz, motherType, histTypeD0DStar, pdgD0);
2785  motherType = 5; FillDStarAndB0Histograms(selectedDStar, primaryVertex, bz, motherType, histTypeDStar);
2786  }
2787 
2788  if(isDesiredCandidate)
2789  {
2790  histType += 1;
2791  motherType = 0; FillD0Histograms(selectedD0, primaryVertex, bz, motherType, histType, pdgD0);
2792  motherType = 1; FillDStarAndB0Histograms(selectedDStar, primaryVertex, bz, motherType, histType);
2793  motherType = 2; FillDStarAndB0Histograms(trackB0, primaryVertex, bz, motherType, histType);
2794  motherType = 3; FillD0Histograms(selectedD0, primaryVertex, bz, motherType, histTypeD0 + 1, pdgD0);
2795  motherType = 4; FillD0Histograms(selectedD0, primaryVertex, bz, motherType, histTypeD0DStar + 1, pdgD0);
2796  motherType = 5; FillDStarAndB0Histograms(selectedDStar, primaryVertex, bz, motherType, histTypeDStar + 1);
2797  }
2798  }
2799  }
2800 
2801  Double_t invmassDelta = DeltaInvMassB0Kpipipi(trackB0);
2802  if(bSameSign)
2803  {
2804  fillthis="invariantMassB0";
2805  fillthis += "_SameSign";
2806  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invariantMassMother);
2807  fillthis="invariantMassB0";
2808  fillthis += ptBinMother + "_SameSign";
2809  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invariantMassMother);
2810  fillthis="invariantMassB0";
2811  fillthis += "_SignSum";
2812  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invariantMassMother,-1);
2813  fillthis="invariantMassB0";
2814  fillthis += ptBinMother + "_SignSum";
2815  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invariantMassMother,-1);
2816  }
2817  if(!bSameSign)
2818  {
2819 
2820  fillthis="deltaInvMassB0";
2821  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invmassDelta);
2822  fillthis="deltaInvMassB0";
2823  fillthis += ptBinMother;
2824  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invmassDelta);
2825 
2826  fillthis="invariantMassB0";
2827  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invariantMassMother);
2828  fillthis="invariantMassB0";
2829  fillthis += ptBinMother;
2830  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invariantMassMother);
2831  fillthis="invariantMassB0";
2832  fillthis += "_SignSum";
2833  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invariantMassMother,1);
2834  fillthis="invariantMassB0";
2835  fillthis += ptBinMother + "_SignSum";
2836  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invariantMassMother,1);
2837  }
2838 
2839  if(trackB0->Pt() > 6.0)
2840  {
2841  ptBinMother = "_ptbin_6_to_inf";
2842  if(bSameSign)
2843  {
2844  fillthis="invariantMassB0";
2845  fillthis += ptBinMother + "_SameSign";
2846  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invariantMassMother);
2847  fillthis="invariantMassB0";
2848  fillthis += ptBinMother + "_SignSum";
2849  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invariantMassMother,-1);
2850  }
2851  if(!bSameSign)
2852  {
2853  fillthis="deltaInvMassB0";
2854  fillthis += ptBinMother;
2855  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invmassDelta);
2856 
2857  fillthis="invariantMassB0";
2858  fillthis += ptBinMother;
2859  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invariantMassMother);
2860  fillthis="invariantMassB0";
2861  fillthis += ptBinMother + "_SignSum";
2862  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invariantMassMother,1);
2863  }
2864 
2865  }
2866 
2867  if(trackB0->Pt() > 3.0)
2868  {
2869  ptBinMother = "_ptbin_3_to_inf";
2870  if(bSameSign)
2871  {
2872  fillthis="invariantMassB0";
2873  fillthis += ptBinMother + "_SameSign";
2874  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invariantMassMother);
2875  fillthis="invariantMassB0";
2876  fillthis += ptBinMother + "_SignSum";
2877  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invariantMassMother,-1);
2878  }
2879  if(!bSameSign)
2880  {
2881  fillthis="deltaInvMassB0";
2882  fillthis += ptBinMother;
2883  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invmassDelta);
2884 
2885  fillthis="invariantMassB0";
2886  fillthis += ptBinMother;
2887  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invariantMassMother);
2888  fillthis="invariantMassB0";
2889  fillthis += ptBinMother + "_SignSum";
2890  ((TH1F*)(fOutputB0MC->FindObject(fillthis)))->Fill(invariantMassMother,1);
2891  }
2892 
2893  }
2894 
2895  delete vertexMother; vertexMother = NULL;
2896  delete trackB0; trackB0 = NULL;
2897  delete vertexDStar; vertexDStar = NULL;
2898  delete trackDStar; trackDStar = NULL;
2899  }
2900  }
2901 
2902  //Restore Object count
2903  //To save space in the table keeping track of all referenced objects,
2904  //we reset the object count to what it was at the beginning of the loop.
2905  TProcessID::SetObjectCount(ObjectNumber);
2906  }
2907  return;
2908 }
2909 //-------------------------------------------------------------------------------------
2910 void AliAnalysisTaskSEB0toDStarPi::FillFinalTrackHistograms(AliAODRecoDecayHF2Prong * selectedB0, Bool_t isDesiredCandidate,TClonesArray * mcTrackArray){
2911 
2912  //In this function we fill histograms with the properties of all the daughters of our selected signal candidate
2913 
2914  AliAODTrack* selectedB0Pion = (AliAODTrack*)selectedB0->GetDaughter(0);
2915  AliAODRecoDecayHF2Prong* selectedDStar = (AliAODRecoDecayHF2Prong*)selectedB0->GetDaughter(1);
2916 
2917  AliAODTrack* selectedDStarPion = (AliAODTrack*)selectedDStar->GetDaughter(0);
2918  AliAODRecoDecayHF2Prong* selectedD0 = (AliAODRecoDecayHF2Prong*)selectedDStar->GetDaughter(1);
2919 
2920  AliAODTrack* selectedD0Pion;
2921  AliAODTrack* selectedD0Kaon;
2922 
2923  if(selectedDStarPion->Charge() == 1) selectedD0Pion = (AliAODTrack*)selectedD0->GetDaughter(0);
2924  if(selectedDStarPion->Charge() == -1) selectedD0Pion = (AliAODTrack*)selectedD0->GetDaughter(1);
2925 
2926  if(selectedDStarPion->Charge() == 1) selectedD0Kaon = (AliAODTrack*)selectedD0->GetDaughter(1);
2927  if(selectedDStarPion->Charge() == -1) selectedD0Kaon = (AliAODTrack*)selectedD0->GetDaughter(0);
2928 
2929  Double_t d0B0pion = TMath::Abs(selectedB0->Getd0Prong(0));
2930  Double_t d0DStarpion = TMath::Abs(selectedDStar->Getd0Prong(0));
2931  Double_t d0D0pion = 0;
2932  Double_t d0D0kaon = 0;
2933 
2934  if(selectedDStarPion->Charge() == 1) d0D0pion = selectedD0->Getd0Prong(0);
2935  if(selectedDStarPion->Charge() == -1) d0D0pion = selectedD0->Getd0Prong(1);
2936 
2937  if(selectedDStarPion->Charge() == 1) d0D0kaon = selectedD0->Getd0Prong(1);
2938  if(selectedDStarPion->Charge() == -1) d0D0kaon = selectedD0->Getd0Prong(0);
2939 
2940  Double_t pt_track = 0;
2941  Double_t momentum_track = 0;
2942  Int_t numberOfITS = 0;
2943  Int_t numberOfTPC = 0;
2944  Int_t daughterType, histType;
2945  Int_t totalNumberOfITS = 0;
2946  Int_t totalNumberOfTPC = 0;
2947  Double_t nSigmaTPC = 0;
2948  Double_t nSigmaTOF = 0;
2949  Double_t nSigmaTPCtotal = 0;
2950  Double_t nSigmaTOFtotal = 0;
2951  Int_t pionPIDnumber = 2;
2952  Int_t kaonPIDnumber = 3;
2953  Int_t TPCok = 0;
2954  Int_t TOFok = 0;
2955 
2956  AliAODPidHF* trackPIDHF = (AliAODPidHF*)fCuts->GetPidHF();
2957 
2958  //fill the D0 pion info
2959  pt_track = selectedD0Pion->Pt();
2960  momentum_track = selectedD0Pion->P();
2961  numberOfITS = selectedD0Pion->GetITSNcls();
2962  numberOfTPC = selectedD0Pion->GetTPCNcls();
2963  totalNumberOfITS += numberOfITS;
2964  totalNumberOfTPC += numberOfTPC;
2965  TPCok = trackPIDHF->GetnSigmaTPC(selectedD0Pion, pionPIDnumber, nSigmaTPC);
2966  TOFok = trackPIDHF->GetnSigmaTOF(selectedD0Pion, pionPIDnumber, nSigmaTOF);
2967  if(TPCok != -1) nSigmaTPCtotal += nSigmaTPC*nSigmaTPC;
2968  if(TOFok != -1) nSigmaTOFtotal += nSigmaTOF*nSigmaTOF;
2969 
2970  Double_t ptB0 = selectedB0->Pt();
2971 
2972  daughterType = 0;
2973  histType = 5;
2974  if(!isDesiredCandidate)
2975  {
2976  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
2977  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
2978  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
2979  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
2980 
2981  for (Int_t j = 0; j < 10; ++j)
2982  {
2983  if(selectedD0Pion->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
2984 
2985  }
2986 
2987  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
2988  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
2989  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
2990  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0D0pion);
2991  ((TH2F*)fDaughterHistogramArray2D[daughterType][4])->Fill(ptB0,pt_track);
2992  }
2993 
2994  if(isDesiredCandidate)
2995  {
2996  histType = 6;
2997  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
2998  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
2999  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
3000  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
3001 
3002  for (Int_t j = 0; j < 10; ++j)
3003  {
3004  if(selectedD0Pion->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
3005 
3006  }
3007 
3008  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
3009  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
3010  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
3011  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0D0pion);
3012  ((TH2F*)fDaughterHistogramArray2D[daughterType][5])->Fill(ptB0,pt_track);
3013  }
3014 
3015  //we save the pdgcode of the used particle and its mother to check PID efficiency
3016  if(fUseMCInfo)
3017  {
3018  Float_t pdgCodeParticle = -1;
3019  Float_t pdgCodeParticleMother = -1;
3020  Int_t mcLabelParticle = -1;
3021  Int_t mcLabelParticleMother = -1;
3022  mcLabelParticle = selectedD0Pion->GetLabel();
3023 
3024  if(mcLabelParticle >= 0){
3025 
3026  AliAODMCParticle *mcTrackParticle = (AliAODMCParticle*)mcTrackArray->At(mcLabelParticle);
3027  pdgCodeParticle = TMath::Abs(mcTrackParticle->GetPdgCode());
3028  ((TH1F*)fDaughterHistogramArray2D[0][2])->Fill(pdgCodeParticle);
3029  mcLabelParticleMother = mcTrackParticle->GetMother();
3030 
3031  if(mcLabelParticleMother >= 0){
3032  AliAODMCParticle *mcTrackParticleMother = (AliAODMCParticle*)mcTrackArray->At(mcLabelParticleMother);
3033  pdgCodeParticleMother = TMath::Abs(mcTrackParticleMother->GetPdgCode());
3034  ((TH1F*)fDaughterHistogramArray2D[0][3])->Fill(pdgCodeParticleMother);
3035  }
3036  }
3037  }
3038 
3039 
3040 
3041  //fill the D0 kaon info
3042  pt_track = selectedD0Kaon->Pt();
3043  momentum_track = selectedD0Kaon->P();
3044  numberOfITS = selectedD0Kaon->GetITSNcls();
3045  numberOfTPC = selectedD0Kaon->GetTPCNcls();
3046  totalNumberOfITS += numberOfITS;
3047  totalNumberOfTPC += numberOfTPC;
3048  TPCok = trackPIDHF->GetnSigmaTPC(selectedD0Kaon, kaonPIDnumber, nSigmaTPC);
3049  TOFok = trackPIDHF->GetnSigmaTOF(selectedD0Kaon, kaonPIDnumber, nSigmaTOF);
3050  if(TPCok != -1) nSigmaTPCtotal += nSigmaTPC*nSigmaTPC;
3051  if(TOFok != -1) nSigmaTOFtotal += nSigmaTOF*nSigmaTOF;
3052 
3053  daughterType = 1;
3054  histType = 5;
3055  if(!isDesiredCandidate)
3056  {
3057  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
3058  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
3059  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
3060  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
3061 
3062  for (Int_t j = 0; j < 10; ++j)
3063  {
3064  if(selectedD0Kaon->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
3065 
3066  }
3067 
3068  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
3069  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
3070  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
3071  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0D0kaon);
3072  ((TH2F*)fDaughterHistogramArray2D[daughterType][4])->Fill(ptB0,pt_track);
3073  }
3074 
3075  if(isDesiredCandidate)
3076  {
3077  histType = 6;
3078  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
3079  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
3080  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
3081  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
3082 
3083  for (Int_t j = 0; j < 10; ++j)
3084  {
3085  if(selectedD0Kaon->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
3086 
3087  }
3088 
3089  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
3090  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
3091  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
3092  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0D0kaon);
3093  ((TH2F*)fDaughterHistogramArray2D[daughterType][5])->Fill(ptB0,pt_track);
3094  }
3095 
3096  //we save the pdgcode of the used particle and its mother to check PID efficiency
3097  if(fUseMCInfo)
3098  {
3099  Float_t pdgCodeParticle = -1;
3100  Float_t pdgCodeParticleMother = -1;
3101  Int_t mcLabelParticle = -1;
3102  Int_t mcLabelParticleMother = -1;
3103  mcLabelParticle = selectedD0Kaon->GetLabel();
3104 
3105  if(mcLabelParticle >= 0){
3106 
3107  AliAODMCParticle *mcTrackParticle = (AliAODMCParticle*)mcTrackArray->At(mcLabelParticle);
3108  pdgCodeParticle = TMath::Abs(mcTrackParticle->GetPdgCode());
3109  ((TH1F*)fDaughterHistogramArray2D[1][2])->Fill(pdgCodeParticle);
3110  mcLabelParticleMother = mcTrackParticle->GetMother();
3111 
3112  if(mcLabelParticleMother >= 0){
3113  AliAODMCParticle *mcTrackParticleMother = (AliAODMCParticle*)mcTrackArray->At(mcLabelParticleMother);
3114  pdgCodeParticleMother = TMath::Abs(mcTrackParticleMother->GetPdgCode());
3115  ((TH1F*)fDaughterHistogramArray2D[1][3])->Fill(pdgCodeParticleMother);
3116  }
3117  }
3118  }
3119 
3120  //fill the DStar pion info
3121  pt_track = selectedDStarPion->Pt();
3122  momentum_track = selectedDStarPion->P();
3123  numberOfITS = selectedDStarPion->GetITSNcls();
3124  numberOfTPC = selectedDStarPion->GetTPCNcls();
3125  totalNumberOfITS += numberOfITS;
3126  totalNumberOfTPC += numberOfTPC;
3127  TPCok = trackPIDHF->GetnSigmaTPC(selectedDStarPion, pionPIDnumber, nSigmaTPC);
3128  TOFok = trackPIDHF->GetnSigmaTOF(selectedDStarPion, pionPIDnumber, nSigmaTOF);
3129  if(TPCok != -1) nSigmaTPCtotal += nSigmaTPC*nSigmaTPC;
3130 
3131  daughterType = 2;
3132  histType = 5;
3133  if(!isDesiredCandidate)
3134  {
3135  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
3136  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
3137  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
3138  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
3139 
3140  for (Int_t j = 0; j < 10; ++j)
3141  {
3142  if(selectedDStarPion->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
3143 
3144  }
3145 
3146  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
3147  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
3148  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
3149  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0DStarpion);
3150  ((TH2F*)fDaughterHistogramArray2D[daughterType][4])->Fill(ptB0,pt_track);
3151  }
3152 
3153  if(isDesiredCandidate)
3154  {
3155  histType = 6;
3156  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
3157  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
3158  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
3159  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
3160 
3161  for (Int_t j = 0; j < 10; ++j)
3162  {
3163  if(selectedDStarPion->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
3164 
3165  }
3166 
3167  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
3168  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
3169  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
3170  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0DStarpion);
3171  ((TH2F*)fDaughterHistogramArray2D[daughterType][5])->Fill(ptB0,pt_track);
3172  }
3173 
3174  //we save the pdgcode of the used particle and its mother to check PID efficiency
3175  if(fUseMCInfo)
3176  {
3177  Float_t pdgCodeParticle = -1;
3178  Float_t pdgCodeParticleMother = -1;
3179  Int_t mcLabelParticle = -1;
3180  Int_t mcLabelParticleMother = -1;
3181  mcLabelParticle = selectedDStarPion->GetLabel();
3182 
3183  if(mcLabelParticle >= 0){
3184 
3185  AliAODMCParticle *mcTrackParticle = (AliAODMCParticle*)mcTrackArray->At(mcLabelParticle);
3186  pdgCodeParticle = TMath::Abs(mcTrackParticle->GetPdgCode());
3187  ((TH1F*)fDaughterHistogramArray2D[2][2])->Fill(pdgCodeParticle);
3188  mcLabelParticleMother = mcTrackParticle->GetMother();
3189 
3190  if(mcLabelParticleMother >= 0){
3191  AliAODMCParticle *mcTrackParticleMother = (AliAODMCParticle*)mcTrackArray->At(mcLabelParticleMother);
3192  pdgCodeParticleMother = TMath::Abs(mcTrackParticleMother->GetPdgCode());
3193  ((TH1F*)fDaughterHistogramArray2D[2][3])->Fill(pdgCodeParticleMother);
3194  }
3195  }
3196  }
3197 
3198  //fill the B0 pion info
3199  pt_track = selectedB0Pion->Pt();
3200  momentum_track = selectedB0Pion->P();
3201  numberOfITS = selectedB0Pion->GetITSNcls();
3202  numberOfTPC = selectedB0Pion->GetTPCNcls();
3203  totalNumberOfITS += numberOfITS;
3204  totalNumberOfTPC += numberOfTPC;
3205  TPCok = trackPIDHF->GetnSigmaTPC(selectedB0Pion, pionPIDnumber, nSigmaTPC);
3206  TOFok = trackPIDHF->GetnSigmaTOF(selectedB0Pion, pionPIDnumber, nSigmaTOF);
3207  if(TPCok != -1) nSigmaTPCtotal += nSigmaTPC*nSigmaTPC;
3208  if(TOFok != -1) nSigmaTOFtotal += nSigmaTOF*nSigmaTOF;
3209 
3210  daughterType = 3;
3211  histType = 5;
3212  if(!isDesiredCandidate)
3213  {
3214  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
3215  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
3216  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
3217  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
3218 
3219  for (Int_t j = 0; j < 10; ++j)
3220  {
3221  if(selectedB0Pion->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
3222 
3223  }
3224 
3225  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
3226  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
3227  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
3228  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0B0pion);
3229  ((TH2F*)fDaughterHistogramArray2D[daughterType][4])->Fill(ptB0,pt_track);
3230  }
3231 
3232  if(isDesiredCandidate)
3233  {
3234  histType = 6;
3235  ((TH1F*)fDaughterHistogramArray[daughterType][histType][0])->Fill(pt_track);
3236  ((TH1F*)fDaughterHistogramArray[daughterType][histType][1])->Fill(momentum_track);
3237  ((TH1F*)fDaughterHistogramArray[daughterType][histType][2])->Fill(numberOfITS);
3238  ((TH1F*)fDaughterHistogramArray[daughterType][histType][3])->Fill(numberOfTPC);
3239 
3240  for (Int_t j = 0; j < 10; ++j)
3241  {
3242  if(selectedB0Pion->HasPointOnITSLayer(j)) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(j);
3243 
3244  }
3245 
3246  if(TPCok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][4])->Fill(nSigmaTPC);
3247  if(TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][5])->Fill(nSigmaTOF);
3248  if(TPCok != -1 && TOFok != -1) ((TH1F*)fDaughterHistogramArray[daughterType][histType][6])->Fill(sqrt(nSigmaTPC*nSigmaTPC + nSigmaTOF*nSigmaTOF));
3249  ((TH1F*)fDaughterHistogramArray[daughterType][histType][8])->Fill(d0B0pion);
3250  ((TH2F*)fDaughterHistogramArray2D[daughterType][5])->Fill(ptB0,pt_track);
3251  }
3252 
3253  //we save the pdgcode of the used particle and its mother to check PID efficiency
3254  if(fUseMCInfo)
3255  {
3256  Float_t pdgCodeParticle = -1;
3257  Float_t pdgCodeParticleMother = -1;
3258  Int_t mcLabelParticle = -1;
3259  Int_t mcLabelParticleMother = -1;
3260  mcLabelParticle = selectedB0Pion->GetLabel();
3261 
3262  if(mcLabelParticle >= 0){
3263 
3264  AliAODMCParticle *mcTrackParticle = (AliAODMCParticle*)mcTrackArray->At(mcLabelParticle);
3265  pdgCodeParticle = TMath::Abs(mcTrackParticle->GetPdgCode());
3266  ((TH1F*)fDaughterHistogramArray2D[3][2])->Fill(pdgCodeParticle);
3267  mcLabelParticleMother = mcTrackParticle->GetMother();
3268 
3269  if(mcLabelParticleMother >= 0){
3270  AliAODMCParticle *mcTrackParticleMother = (AliAODMCParticle*)mcTrackArray->At(mcLabelParticleMother);
3271  pdgCodeParticleMother = TMath::Abs(mcTrackParticleMother->GetPdgCode());
3272  ((TH1F*)fDaughterHistogramArray2D[3][3])->Fill(pdgCodeParticleMother);
3273  }
3274  }
3275  }
3276 
3277  if(!isDesiredCandidate)
3278  {
3279  ((TH1F*)(fOutputB0MC->FindObject("totalITSBackground")))->Fill(totalNumberOfITS);
3280  ((TH1F*)(fOutputB0MC->FindObject("totalTPCBackground")))->Fill(totalNumberOfTPC);
3281  ((TH1F*)(fOutputB0MC->FindObject("totalSigmaPIDBackground")))->Fill(sqrt(nSigmaTPCtotal + nSigmaTOFtotal));
3282  }
3283  if(isDesiredCandidate)
3284  {
3285  ((TH1F*)(fOutputB0MC->FindObject("totalITSSignal")))->Fill(totalNumberOfITS);
3286  ((TH1F*)(fOutputB0MC->FindObject("totalTPCSignal")))->Fill(totalNumberOfTPC);
3287  ((TH1F*)(fOutputB0MC->FindObject("totalSigmaPIDSignal")))->Fill(sqrt(nSigmaTPCtotal + nSigmaTOFtotal));
3288  }
3289 
3290 
3291  return;
3292 }
3293 //-------------------------------------------------------------------------------------
3295 {
3299  Int_t chargeDStar = DStar->Charge();
3300 
3301  Double_t e[3];
3302  UInt_t prongs[2];
3303  if(chargeDStar==1){
3304  e[0]=((AliAODRecoDecayHF2Prong*)DStar->GetDaughter(1))->EProng(0,211);
3305  e[1]=((AliAODRecoDecayHF2Prong*)DStar->GetDaughter(1))->EProng(1,321);
3306  prongs[0] = 211;
3307  prongs[1] = 321;
3308  }
3309  else if (chargeDStar==-1)
3310  {
3311  e[0]=((AliAODRecoDecayHF2Prong*)DStar->GetDaughter(1))->EProng(1,211);
3312  e[1]=((AliAODRecoDecayHF2Prong*)DStar->GetDaughter(1))->EProng(0,321);
3313  prongs[1] = 211;
3314  prongs[0] = 321;
3315  }
3316  else
3317  {
3318  std::cout << "Wrong charge DStar." << std::endl;
3319  return 0;
3320  }
3321  e[2]=DStar->EProng(0,211);
3322 
3323  Double_t esum = e[0]+e[1]+e[2];
3324  Double_t invMassDStar = TMath::Sqrt(esum*esum-DStar->P()*DStar->P());
3325  Double_t invMassD0 = ((AliAODRecoDecayHF2Prong*)DStar->GetDaughter(1))->InvMass(2,prongs);
3326 
3327  return invMassDStar - invMassD0;
3328 }
3329 //-------------------------------------------------------------------------------------
3331 {
3335 
3336  AliAODRecoDecayHF2Prong * DStar = (AliAODRecoDecayHF2Prong*)B0->GetDaughter(1);
3337  Int_t chargeDStar = DStar->Charge();
3338 
3339  Double_t e[4];
3340  UInt_t prongs[2];
3341  if(chargeDStar==1){
3342  e[1]=((AliAODRecoDecayHF2Prong*)DStar->GetDaughter(1))->EProng(0,211);
3343  e[2]=((AliAODRecoDecayHF2Prong*)DStar->GetDaughter(1))->EProng(1,321);
3344  prongs[0] = 211;
3345  prongs[1] = 321;
3346  }
3347  else if (chargeDStar==-1)
3348  {
3349  e[1]=((AliAODRecoDecayHF2Prong*)DStar->GetDaughter(1))->EProng(1,211);
3350  e[2]=((AliAODRecoDecayHF2Prong*)DStar->GetDaughter(1))->EProng(0,321);
3351  prongs[1] = 211;
3352  prongs[0] = 321;
3353  }
3354  else
3355  {
3356  std::cout << "Wrong charge DStar." << std::endl;
3357  return 0;
3358  }
3359  e[0]=DStar->EProng(0,211);
3360  e[3]=B0->EProng(0,211);
3361 
3362  Double_t esum = e[0]+e[1]+e[2]+e[3];
3363  Double_t invMassB0 = TMath::Sqrt(esum*esum-B0->P()*B0->P());
3364 
3365  Double_t invMassD0 = ((AliAODRecoDecayHF2Prong*)DStar->GetDaughter(1))->InvMass(2,prongs);
3366 
3367  return invMassB0 - invMassD0;
3368 }
3369 //-------------------------------------------------------------------------------------
3370 void AliAnalysisTaskSEB0toDStarPi::FillD0Histograms(AliAODRecoDecayHF2Prong * selectedMother, AliAODVertex *primaryVertex, Double_t bz, Int_t motherType, Int_t histType, Int_t pdgCodeMother){
3371 
3372  if(histType<0) return;
3373 
3374  //In this function we fill the histograms of the reconstructed mothers
3375  Double_t ptMother = 0.0;
3376  Double_t momentumMother = 0.0;
3377  Double_t etaMother = 0.0;
3378  Double_t phiMother = 0.0;
3379  Double_t d0Mother = 0.0;
3380  Double_t d0firstTrack = 0.0;
3381  Double_t d0secondTrack = 0.0;
3382  Double_t pointingAngle = 0.0;
3383  Double_t impactProduct = 0.0;
3384  Double_t impactProductXY = 0.0;
3385  Double_t invariantMassMother = 0.0;
3386  Double_t invmassDelta = 0.0;
3387  Double_t dcaMother = 0.0;
3388  AliAODVertex * vertexMother = 0x0;
3389  AliAODVertex * vertexDaughter = 0x0;
3390  Double_t vertexDistance = 0.0;
3391  Double_t decayLengthDaughter = 0.0;
3392  Double_t decayTime = 0.0;
3393  Double_t angleMotherFirstDaughter = 0.0;
3394  Double_t angleMotherSecondDaughter = 0.0;
3395  Double_t topomaticFirstDaughter = 0.0;
3396  Double_t topomaticSecondDaughter = 0.0;
3397  Double_t ptFirstDaughter = 0.0;
3398  Double_t ptSecondDaughter = 0.0;
3399  UInt_t prongs[2], prongs2[2];
3400  Double_t angleBetweenBothDaughters = 0;
3401  Double_t cosThetaStar = 0;
3402  Double_t normDecayLength = 0;
3403  Double_t pdgMassMother=TDatabasePDG::Instance()->GetParticle(421)->Mass();
3404 
3405 
3406  prongs[0] = 211; prongs[1] = 321;
3407  prongs2[0] = 321; prongs2[1] = 211;
3408  AliAODTrack * firstDaughter = (AliAODTrack*)selectedMother->GetDaughter(0);
3409  AliAODTrack * secondDaughter = (AliAODTrack*)selectedMother->GetDaughter(1);
3410  vertexMother = selectedMother->GetSecondaryVtx();
3411  ptFirstDaughter = firstDaughter->Pt();
3412  ptSecondDaughter = secondDaughter->Pt();
3413 
3414  //Topomatic
3415  Double_t dd0pr1=0.;
3416  Double_t dd0pr2=0.;
3417  Double_t dd0max=0.;
3418  Double_t dd0min=0.;
3419  for(Int_t ipr=0; ipr<2; ipr++)
3420  {
3421  Double_t diffIP, errdiffIP;
3422  selectedMother->Getd0MeasMinusExpProng(ipr,bz,diffIP,errdiffIP);
3423  Double_t normdd0=0.;
3424  if(errdiffIP>0.) normdd0=diffIP/errdiffIP;
3425  if(ipr==0) dd0pr1=normdd0;
3426  if(ipr==1) dd0pr2=normdd0;
3427 
3428  }
3429  if(TMath::Abs(dd0pr1)>TMath::Abs(dd0pr2)) {dd0max=dd0pr1; dd0min=dd0pr2;}
3430  else {dd0max=dd0pr2; dd0min=dd0pr1;}
3431 
3432  AliExternalTrackParam motherTrack;
3433  motherTrack.CopyFromVTrack(selectedMother);
3434  Double_t d0z0[2],covd0z0[3],d0[2];
3435  motherTrack.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
3436  d0[0] = d0z0[0];
3437 
3438  ptMother = selectedMother->Pt();
3439  momentumMother = selectedMother->P();
3440  etaMother = selectedMother->Eta();
3441  phiMother = selectedMother->Phi();
3442 
3443  d0Mother = TMath::Abs(d0[0]);
3444  d0firstTrack = TMath::Abs(selectedMother->Getd0Prong(0));
3445  d0secondTrack = TMath::Abs(selectedMother->Getd0Prong(1));
3446  pointingAngle = selectedMother->CosPointingAngle();
3447  impactProduct = selectedMother->Prodd0d0();
3448  impactProductXY = TMath::Abs(selectedMother->ImpParXY());
3449  invariantMassMother = selectedMother->InvMass(2,prongs);
3450  if(pdgCodeMother == -421) invariantMassMother = selectedMother->InvMass(2,prongs2);
3451 
3452  dcaMother = selectedMother->GetDCA();
3453  vertexDistance = vertexMother->DistanceToVertex(primaryVertex);
3454  angleMotherFirstDaughter = (selectedMother->Px() * firstDaughter->Px() + selectedMother->Py() * firstDaughter->Py() + selectedMother->Pz() * firstDaughter->Pz()) /(selectedMother->P() * firstDaughter->P());
3455  angleMotherSecondDaughter = (selectedMother->Px() * secondDaughter->Px() + selectedMother->Py() * secondDaughter->Py() + selectedMother->Pz() * secondDaughter->Pz()) /(selectedMother->P() * secondDaughter->P());
3456  cosThetaStar = selectedMother->CosThetaStar(0,421,211,321);
3457  angleBetweenBothDaughters = (firstDaughter->Px() * secondDaughter->Px() + firstDaughter->Py() * secondDaughter->Py() + firstDaughter->Pz() * secondDaughter->Pz()) /(firstDaughter->P() * secondDaughter->P());
3458  normDecayLength = selectedMother->NormalizedDecayLength();
3459 
3460  Double_t pseudoProperDecayLength = ((vertexMother->GetX() - primaryVertex->GetX()) * selectedMother->Px() / TMath::Abs(selectedMother->Pt())) + ((vertexMother->GetY() - primaryVertex->GetY()) * selectedMother->Py() / TMath::Abs(selectedMother->Pt()));
3461  Double_t pseudoProperDecayTime = pseudoProperDecayLength * pdgMassMother/ptMother;
3462  decayTime = vertexDistance / (299792458 * TMath::Sqrt(1/((pdgMassMother*pdgMassMother/(momentumMother*momentumMother)) + 1)));
3463 
3464  Double_t phi = selectedMother->Phi();
3465  Double_t theta = selectedMother->Theta();
3466  Double_t covMatrix[21];
3467  selectedMother->GetCovarianceXYZPxPyPz(covMatrix);
3468 
3469  Double_t cp = TMath::Cos(phi);
3470  Double_t sp = TMath::Sin(phi);
3471  Double_t ct = TMath::Cos(theta);
3472  Double_t st = TMath::Sin(theta);
3473 
3474  Double_t errorMomentum = covMatrix[9]*cp*cp*ct*ct // GetCovPxPx
3475  +covMatrix[13]*2.*cp*sp*ct*ct // GetCovPxPy
3476  +covMatrix[18]*2.*cp*ct*st // GetCovPxPz
3477  +covMatrix[14]*sp*sp*ct*ct // GetCovPyPy
3478  +covMatrix[19]*2.*sp*ct*st // GetCovPyPz
3479  +covMatrix[20]*st*st; // GetCovPzPz
3480  Double_t normalizedDecayTime = selectedMother->NormalizedDecayLength() / (299792458 * TMath::Sqrt(1/((pdgMassMother*pdgMassMother*errorMomentum*errorMomentum/(momentumMother*momentumMother)) + 1)));
3481 
3482  Double_t eKaon = selectedMother->EProng(1,321);
3483  Double_t invMassKaon = TMath::Sqrt(eKaon*eKaon-secondDaughter->P()*secondDaughter->P());
3484  Double_t invMassD0 = selectedMother->InvMassD0();
3485  invmassDelta = invMassD0 - invMassKaon;
3486 
3487  Double_t vertexMotherX = vertexMother->GetX();
3488  Double_t vertexMotherY = vertexMother->GetY();
3489  Double_t vertexMotherZ = vertexMother->GetZ();
3490 
3491  ((TH1F*)fMotherHistogramArray[motherType][histType][0])->Fill(ptMother);
3492  ((TH1F*)fMotherHistogramArray[motherType][histType][1])->Fill(ptFirstDaughter);
3493  ((TH1F*)fMotherHistogramArray[motherType][histType][2])->Fill(ptSecondDaughter);
3494  ((TH1F*)fMotherHistogramArray[motherType][histType][3])->Fill(etaMother);
3495  ((TH1F*)fMotherHistogramArray[motherType][histType][4])->Fill(phiMother);
3496  ((TH1F*)fMotherHistogramArray[motherType][histType][5])->Fill(d0Mother);
3497  ((TH1F*)fMotherHistogramArray[motherType][histType][6])->Fill(d0firstTrack);
3498  ((TH1F*)fMotherHistogramArray[motherType][histType][7])->Fill(d0secondTrack);
3499  ((TH1F*)fMotherHistogramArray[motherType][histType][8])->Fill(pointingAngle);
3500  ((TH1F*)fMotherHistogramArray[motherType][histType][9])->Fill(impactProduct);
3501  ((TH1F*)fMotherHistogramArray[motherType][histType][10])->Fill(impactProductXY);
3502  ((TH1F*)fMotherHistogramArray[motherType][histType][11])->Fill(invariantMassMother);
3503  ((TH1F*)fMotherHistogramArray[motherType][histType][12])->Fill(invmassDelta);
3504  ((TH1F*)fMotherHistogramArray[motherType][histType][13])->Fill(dcaMother);
3505  ((TH1F*)fMotherHistogramArray[motherType][histType][14])->Fill(vertexDistance);
3506  ((TH1F*)fMotherHistogramArray[motherType][histType][15])->Fill(normDecayLength);
3507  ((TH1F*)fMotherHistogramArray[motherType][histType][16])->Fill(pseudoProperDecayTime);
3508  ((TH1F*)fMotherHistogramArray[motherType][histType][17])->Fill(decayTime);
3509  ((TH1F*)fMotherHistogramArray[motherType][histType][18])->Fill(normalizedDecayTime);
3510  ((TH1F*)fMotherHistogramArray[motherType][histType][19])->Fill(angleMotherFirstDaughter);
3511  ((TH1F*)fMotherHistogramArray[motherType][histType][20])->Fill(angleMotherSecondDaughter);
3512  ((TH1F*)fMotherHistogramArray[motherType][histType][21])->Fill(angleBetweenBothDaughters);
3513  ((TH1F*)fMotherHistogramArray[motherType][histType][22])->Fill(cosThetaStar);
3514 
3515  ((TH1F*)fMotherHistogramArray[motherType][histType][23])->Fill(vertexMotherX);
3516  ((TH1F*)fMotherHistogramArray[motherType][histType][24])->Fill(vertexMotherY);
3517  ((TH1F*)fMotherHistogramArray[motherType][histType][25])->Fill(vertexMotherZ);
3518 
3519  ((TH1F*)fMotherHistogramArray[motherType][histType][36])->Fill(TMath::Abs(dd0pr1));
3520  ((TH1F*)fMotherHistogramArray[motherType][histType][37])->Fill(TMath::Abs(dd0pr2));
3521  ((TH1F*)fMotherHistogramArray[motherType][histType][38])->Fill(TMath::Abs(dd0max));
3522  ((TH1F*)fMotherHistogramArray[motherType][histType][39])->Fill(TMath::Abs(dd0min));
3523 
3524 
3525 
3526  //we fill the 2D histograms
3527  Int_t nFirst = 0;
3528  Int_t nSecond = 1;
3529  Int_t nVariables = 8;
3530  Int_t nHistograms = nVariables * (nVariables - 1) / 2;
3531  for (Int_t k = 0; k < nHistograms; ++k)
3532  {
3533  Double_t firstVariable = 0.0;
3534  Double_t secondVariable = 0.0;
3535 
3536  if(nFirst==0) firstVariable = d0firstTrack;
3537  if(nFirst==1) firstVariable = d0secondTrack;
3538  if(nFirst==2) firstVariable = d0Mother;
3539  if(nFirst==3) firstVariable = pointingAngle;
3540  if(nFirst==4) firstVariable = impactProduct;
3541  if(nFirst==5) firstVariable = impactProductXY;
3542  if(nFirst==6) firstVariable = vertexDistance;
3543  if(nFirst==7) firstVariable = normDecayLength;
3544  if(nFirst==8) firstVariable = pseudoProperDecayTime;
3545 
3546  if(nSecond==0) secondVariable = d0firstTrack;
3547  if(nSecond==1) secondVariable = d0secondTrack;
3548  if(nSecond==2) secondVariable = d0Mother;
3549  if(nSecond==3) secondVariable = pointingAngle;
3550  if(nSecond==4) secondVariable = impactProduct;
3551  if(nSecond==5) secondVariable = impactProductXY;
3552  if(nSecond==6) secondVariable = vertexDistance;
3553  if(nSecond==7) secondVariable = normDecayLength;
3554  if(nSecond==8) secondVariable = pseudoProperDecayTime;
3555 
3556  ((TH2F*)fMotherHistogramArray2D[motherType][histType][k])->Fill(firstVariable,secondVariable);
3557 
3558  nSecond++;
3559  if(nSecond>nVariables)
3560  {
3561  nFirst++;
3562  nSecond = nFirst + 1;
3563  }
3564  }
3565 
3566  return;
3567 }
3568 //-------------------------------------------------------------------------------------
3569 void AliAnalysisTaskSEB0toDStarPi::FillDStarAndB0Histograms(AliAODRecoDecayHF2Prong * selectedMother, AliAODVertex *primaryVertex, Double_t bz, Int_t motherType, Int_t histType){
3570 
3571  if(histType<0) return;
3572 
3573  //In this function we fill the histograms of the reconstructed mothers
3574  Double_t ptMother = 0.0;
3575  Double_t momentumMother = 0.0;
3576  Double_t etaMother = 0.0;
3577  Double_t phiMother = 0.0;
3578  Double_t d0Mother = 0.0;
3579  Double_t d0firstTrack = 0.0;
3580  Double_t d0secondTrack = 0.0;
3581  Double_t pointingAngle = 0.0;
3582  Double_t impactProduct = 0.0;
3583  Double_t impactProductXY = 0.0;
3584  Double_t invariantMassMother = 0.0;
3585  Double_t invmassDelta = 0.0;
3586  Double_t dcaMother = 0.0;
3587  AliAODVertex * vertexMother = 0x0;
3588  AliAODVertex * vertexDaughter = 0x0;
3589  Double_t vertexDistance = 0.0;
3590  Double_t normDecayLength = 0.0;
3591  Double_t decayTime = 0.0;
3592  Double_t angleMotherFirstDaughter = 0.0;
3593  Double_t angleMotherSecondDaughter = 0.0;
3594  Double_t topomaticFirstDaughter = 0.0;
3595  Double_t topomaticSecondDaughter = 0.0;
3596  Double_t ptFirstDaughter = 0.0;
3597  Double_t ptSecondDaughter = 0.0;
3598  UInt_t prongs[2];
3599  Double_t cosThetaStar = 0;
3600  Double_t angleBetweenBothDaughters = 0;
3601 
3602  Double_t pdgMassMother = 0;
3603 
3604  if(motherType==1 || motherType==5){ //DStar
3605  prongs[0] = 211; prongs[1] = 421;
3606  invmassDelta = DeltaInvMassDStarKpipi(selectedMother);
3607  AliAODTrack * firstDaughter = (AliAODTrack*)selectedMother->GetDaughter(0);
3608  AliAODRecoDecayHF2Prong * secondDaughter = (AliAODRecoDecayHF2Prong*)selectedMother->GetDaughter(1);
3609 
3610  ptFirstDaughter = firstDaughter->Pt();
3611  ptSecondDaughter = secondDaughter->Pt();
3612  vertexMother = selectedMother->GetSecondaryVtx();
3613  vertexDaughter = secondDaughter->GetSecondaryVtx();
3614  angleMotherFirstDaughter = (selectedMother->Px() * firstDaughter->Px() + selectedMother->Py() * firstDaughter->Py() + selectedMother->Pz() * firstDaughter->Pz()) /(selectedMother->P() * firstDaughter->P());
3615  angleMotherSecondDaughter = (selectedMother->Px() * secondDaughter->Px() + selectedMother->Py() * secondDaughter->Py() + selectedMother->Pz() * secondDaughter->Pz()) /(selectedMother->P() * secondDaughter->P());
3616  angleBetweenBothDaughters = (firstDaughter->Px() * secondDaughter->Px() + firstDaughter->Py() * secondDaughter->Py() + firstDaughter->Pz() * secondDaughter->Pz()) /(firstDaughter->P() * secondDaughter->P());
3617  cosThetaStar = selectedMother->CosThetaStar(0,413,211,421);
3618  pdgMassMother = TDatabasePDG::Instance()->GetParticle(413)->Mass();
3619  }
3620  if(motherType==2){ //B0
3621  prongs[0] = 211; prongs[1] = 413;
3622  invmassDelta = DeltaInvMassB0Kpipipi(selectedMother);
3623  AliAODTrack * firstDaughter = (AliAODTrack*)selectedMother->GetDaughter(0);
3624  AliAODRecoDecayHF2Prong * secondDaughter = (AliAODRecoDecayHF2Prong*)selectedMother->GetDaughter(1);
3625 
3626  ptFirstDaughter = firstDaughter->Pt();
3627  ptSecondDaughter = secondDaughter->Pt();
3628  vertexMother = selectedMother->GetSecondaryVtx();
3629  vertexDaughter = secondDaughter->GetSecondaryVtx();
3630  angleMotherFirstDaughter = (selectedMother->Px() * firstDaughter->Px() + selectedMother->Py() * firstDaughter->Py() + selectedMother->Pz() * firstDaughter->Pz()) /(selectedMother->P() * firstDaughter->P());
3631  angleMotherSecondDaughter = (selectedMother->Px() * secondDaughter->Px() + selectedMother->Py() * secondDaughter->Py() + selectedMother->Pz() * secondDaughter->Pz()) /(selectedMother->P() * secondDaughter->P());
3632  angleBetweenBothDaughters = (firstDaughter->Px() * secondDaughter->Px() + firstDaughter->Py() * secondDaughter->Py() + firstDaughter->Pz() * secondDaughter->Pz()) /(firstDaughter->P() * secondDaughter->P());
3633  cosThetaStar = selectedMother->CosThetaStar(0,511,211,413);
3634  pdgMassMother = TDatabasePDG::Instance()->GetParticle(511)->Mass();
3635  }
3636 
3637 
3638  //Topomatic
3639  Double_t dd0pr1=0.;
3640  Double_t dd0pr2=0.;
3641  Double_t dd0max=0.;
3642  Double_t dd0min=0.;
3643  for(Int_t ipr=0; ipr<2; ipr++)
3644  {
3645  Double_t diffIP, errdiffIP;
3646  selectedMother->Getd0MeasMinusExpProng(ipr,bz,diffIP,errdiffIP);
3647  Double_t normdd0=0.;
3648  if(errdiffIP>0.) normdd0=diffIP/errdiffIP;
3649  if(ipr==0) dd0pr1=normdd0;
3650  if(ipr==1) dd0pr2=normdd0;
3651 
3652  // else if(TMath::Abs(normdd0)>TMath::Abs(dd0max)) dd0max=normdd0;
3653  }
3654  if(TMath::Abs(dd0pr1)>TMath::Abs(dd0pr2)) {dd0max=dd0pr1; dd0min=dd0pr2;}
3655  else {dd0max=dd0pr2; dd0min=dd0pr1;}
3656 
3657 
3658 
3659  AliExternalTrackParam motherTrack;
3660  motherTrack.CopyFromVTrack(selectedMother);
3661  Double_t d0z0[2],covd0z0[3],d0[2];
3662  motherTrack.PropagateToDCA(primaryVertex,bz,100.,d0z0,covd0z0);
3663  d0[0] = d0z0[0];
3664 
3665  ptMother = selectedMother->Pt();
3666  momentumMother = selectedMother->P();
3667  etaMother = selectedMother->Eta();
3668  phiMother = selectedMother->Phi();
3669  d0Mother = TMath::Abs(d0[0]);
3670 
3671  pointingAngle = selectedMother->CosPointingAngle();
3672  impactProduct = selectedMother->Prodd0d0();
3673  impactProductXY = TMath::Abs(selectedMother->ImpParXY());
3674  invariantMassMother = selectedMother->InvMass(2,prongs);
3675  dcaMother = selectedMother->GetDCA();
3676  vertexDistance = vertexMother->DistanceToVertex(primaryVertex);
3677  d0firstTrack = TMath::Abs(selectedMother->Getd0Prong(0));
3678  d0secondTrack = TMath::Abs(selectedMother->Getd0Prong(1));
3679  normDecayLength = selectedMother->NormalizedDecayLength();
3680 
3681  Double_t pseudoProperDecayLength = ((vertexMother->GetX() - primaryVertex->GetX()) * selectedMother->Px() / TMath::Abs(selectedMother->Pt())) + ((vertexMother->GetY() - primaryVertex->GetY()) * selectedMother->Py() / TMath::Abs(selectedMother->Pt()));
3682  Double_t pseudoProperDecayTime = pseudoProperDecayLength * pdgMassMother/ptMother;
3683  decayTime = vertexDistance / (299792458 * TMath::Sqrt(1/((pdgMassMother*pdgMassMother/(momentumMother*momentumMother)) + 1)));
3684 
3685  Double_t phi = selectedMother->Phi();
3686  Double_t theta = selectedMother->Theta();
3687  Double_t covMatrix[21];
3688  selectedMother->GetCovarianceXYZPxPyPz(covMatrix);
3689 
3690  Double_t cp = TMath::Cos(phi);
3691  Double_t sp = TMath::Sin(phi);
3692  Double_t ct = TMath::Cos(theta);
3693  Double_t st = TMath::Sin(theta);
3694 
3695  Double_t errorMomentum = covMatrix[9]*cp*cp*ct*ct // GetCovPxPx
3696  +covMatrix[13]*2.*cp*sp*ct*ct // GetCovPxPy
3697  +covMatrix[18]*2.*cp*ct*st // GetCovPxPz
3698  +covMatrix[14]*sp*sp*ct*ct // GetCovPyPy
3699  +covMatrix[19]*2.*sp*ct*st // GetCovPyPz
3700  +covMatrix[20]*st*st; // GetCovPzPz
3701  Double_t normalizedDecayTime = selectedMother->NormalizedDecayLength() / (299792458 * TMath::Sqrt(1/((pdgMassMother*pdgMassMother*errorMomentum*errorMomentum/(momentumMother*momentumMother)) + 1)));
3702 
3703  Double_t vertexMotherX = vertexMother->GetX();
3704  Double_t vertexMotherY = vertexMother->GetY();
3705  Double_t vertexMotherZ = vertexMother->GetZ();
3706 
3707  ((TH1F*)fMotherHistogramArray[motherType][histType][0])->Fill(ptMother);
3708  ((TH1F*)fMotherHistogramArray[motherType][histType][1])->Fill(ptFirstDaughter);
3709  ((TH1F*)fMotherHistogramArray[motherType][histType][2])->Fill(ptSecondDaughter);
3710  ((TH1F*)fMotherHistogramArray[motherType][histType][3])->Fill(etaMother);
3711  ((TH1F*)fMotherHistogramArray[motherType][histType][4])->Fill(phiMother);
3712  ((TH1F*)fMotherHistogramArray[motherType][histType][5])->Fill(d0Mother);
3713  ((TH1F*)fMotherHistogramArray[motherType][histType][6])->Fill(d0firstTrack);
3714  ((TH1F*)fMotherHistogramArray[motherType][histType][7])->Fill(d0secondTrack);
3715  ((TH1F*)fMotherHistogramArray[motherType][histType][8])->Fill(pointingAngle);
3716  ((TH1F*)fMotherHistogramArray[motherType][histType][9])->Fill(impactProduct);
3717  ((TH1F*)fMotherHistogramArray[motherType][histType][10])->Fill(impactProductXY);
3718  ((TH1F*)fMotherHistogramArray[motherType][histType][11])->Fill(invariantMassMother);
3719  ((TH1F*)fMotherHistogramArray[motherType][histType][12])->Fill(invmassDelta);
3720  ((TH1F*)fMotherHistogramArray[motherType][histType][13])->Fill(dcaMother);
3721  ((TH1F*)fMotherHistogramArray[motherType][histType][14])->Fill(vertexDistance);
3722  ((TH1F*)fMotherHistogramArray[motherType][histType][15])->Fill(normDecayLength);
3723  ((TH1F*)fMotherHistogramArray[motherType][histType][16])->Fill(pseudoProperDecayTime);
3724  ((TH1F*)fMotherHistogramArray[motherType][histType][17])->Fill(decayTime);
3725  ((TH1F*)fMotherHistogramArray[motherType][histType][18])->Fill(normalizedDecayTime);
3726  ((TH1F*)fMotherHistogramArray[motherType][histType][19])->Fill(angleMotherFirstDaughter);
3727  ((TH1F*)fMotherHistogramArray[motherType][histType][20])->Fill(angleMotherSecondDaughter);
3728  ((TH1F*)fMotherHistogramArray[motherType][histType][21])->Fill(angleBetweenBothDaughters);
3729  ((TH1F*)fMotherHistogramArray[motherType][histType][22])->Fill(cosThetaStar);
3730 
3731  ((TH1F*)fMotherHistogramArray[motherType][histType][23])->Fill(vertexMotherX);
3732  ((TH1F*)fMotherHistogramArray[motherType][histType][24])->Fill(vertexMotherY);
3733  ((TH1F*)fMotherHistogramArray[motherType][histType][25])->Fill(vertexMotherZ);
3734 
3735  ((TH1F*)fMotherHistogramArray[motherType][histType][36])->Fill(TMath::Abs(dd0pr1));
3736  ((TH1F*)fMotherHistogramArray[motherType][histType][37])->Fill(TMath::Abs(dd0pr2));
3737  ((TH1F*)fMotherHistogramArray[motherType][histType][38])->Fill(TMath::Abs(dd0max));
3738  ((TH1F*)fMotherHistogramArray[motherType][histType][39])->Fill(TMath::Abs(dd0min));
3739 
3740  //we fill the 2D histograms
3741  Int_t nFirst = 0;
3742  Int_t nSecond = 1;
3743  Int_t nVariables = 8;
3744  Int_t nHistograms = nVariables * (nVariables - 1) / 2;
3745  for (Int_t k = 0; k < nHistograms; ++k)
3746  {
3747  Double_t firstVariable = 0.0;
3748  Double_t secondVariable = 0.0;
3749 
3750  if(nFirst==0) firstVariable = d0firstTrack;
3751  if(nFirst==1) firstVariable = d0secondTrack;
3752  if(nFirst==2) firstVariable = d0Mother;
3753  if(nFirst==3) firstVariable = pointingAngle;
3754  if(nFirst==4) firstVariable = impactProduct;
3755  if(nFirst==5) firstVariable = impactProductXY;
3756  if(nFirst==6) firstVariable = vertexDistance;
3757  if(nFirst==7) firstVariable = normDecayLength;
3758  if(nFirst==8) firstVariable = pseudoProperDecayTime;
3759 
3760  if(nSecond==0) secondVariable = d0firstTrack;
3761  if(nSecond==1) secondVariable = d0secondTrack;
3762  if(nSecond==2) secondVariable = d0Mother;
3763  if(nSecond==3) secondVariable = pointingAngle;
3764  if(nSecond==4) secondVariable = impactProduct;
3765  if(nSecond==5) secondVariable = impactProductXY;
3766  if(nSecond==6) secondVariable = vertexDistance;
3767  if(nSecond==7) secondVariable = normDecayLength;
3768  if(nSecond==8) secondVariable = pseudoProperDecayTime;
3769 
3770  ((TH2F*)fMotherHistogramArray2D[motherType][histType][k])->Fill(firstVariable,secondVariable);
3771 
3772  nSecond++;
3773  if(nSecond>nVariables)
3774  {
3775  nFirst++;
3776  nSecond = nFirst + 1;
3777  }
3778  }
3779 
3780  if(motherType==1 || motherType==5){
3781  motherType = motherType -1;
3782 
3783  AliAODRecoDecay* secondDaughter = (AliAODRecoDecay*)selectedMother->GetDaughter(1);
3784  AliAODTrack * pionD0 = (AliAODTrack*)selectedMother->GetDaughter(0);
3785  AliAODTrack * kaonD0 = (AliAODTrack*)selectedMother->GetDaughter(1);
3786 
3787  AliAODVertex * vertexDStar = vertexMother;
3788  AliAODVertex * vertexD0 = secondDaughter->GetSecondaryVtx();
3789  pdgMassMother = TDatabasePDG::Instance()->GetParticle(421)->Mass();
3790 
3791  AliExternalTrackParam pionD0Track;
3792  AliExternalTrackParam kaonD0Track;
3793 
3794  Double_t d0z0[2],covd0z0[3],d0[2];
3795 
3796  pionD0Track.CopyFromVTrack(pionD0);
3797  pionD0Track.PropagateToDCA(vertexDStar,bz,100.,d0z0,covd0z0);
3798  d0[0] = d0z0[0];
3799 
3800  kaonD0Track.CopyFromVTrack(kaonD0);
3801  kaonD0Track.PropagateToDCA(vertexDStar,bz,100.,d0z0,covd0z0);
3802  d0[1] = d0z0[0];
3803 
3804  AliExternalTrackParam D0Track;
3805  D0Track.CopyFromVTrack(secondDaughter);
3806  Double_t d0z0D0[2],covd0z0D0[3],d0D0;
3807  motherTrack.PropagateToDCA(vertexDStar,bz,100.,d0z0D0,covd0z0D0);
3808  d0D0 = d0z0D0[0];
3809 
3810  Double_t impactProductToDStar = d0[0]*d0[1];
3811  Double_t impactProductXYToDStar = secondDaughter->ImpParXY(vertexDStar);
3812 
3813  Double_t momentumMother = secondDaughter->P();
3814  Double_t pointingAngleToDStar = secondDaughter->CosPointingAngle(vertexDStar);
3815  Double_t d0FirstDaughterToDStar = TMath::Abs(d0[0]);
3816  Double_t d0SecondDaughterToDStar = TMath::Abs(d0[1]);
3817  Double_t normDecayLengthToDStar = secondDaughter->NormalizedDecayLength(vertexDStar);
3818 
3819  Double_t pseudoProperDecayLength = ((vertexD0->GetX() - vertexDStar->GetX()) * secondDaughter->Px() / TMath::Abs(secondDaughter->Pt())) + ((vertexD0->GetY() - vertexDStar->GetY()) * secondDaughter->Py() / TMath::Abs(secondDaughter->Pt()));
3820  Double_t pseudoProperDecayTimeToDStar = pseudoProperDecayLength * pdgMassMother/ptMother;
3821  Double_t DecayTimeToDStar = vertexDistance / (299792458 * TMath::Sqrt(1/((pdgMassMother*pdgMassMother/(momentumMother*momentumMother)) + 1)));
3822 
3823  Double_t phi = secondDaughter->Phi();
3824  Double_t theta = secondDaughter->Theta();
3825  Double_t covMatrix[21];
3826  secondDaughter->GetCovarianceXYZPxPyPz(covMatrix);
3827 
3828  Double_t cp = TMath::Cos(phi);
3829  Double_t sp = TMath::Sin(phi);
3830  Double_t ct = TMath::Cos(theta);
3831  Double_t st = TMath::Sin(theta);
3832 
3833  Double_t errorMomentum = covMatrix[9]*cp*cp*ct*ct // GetCovPxPx
3834  +covMatrix[13]*2.*cp*sp*ct*ct // GetCovPxPy
3835  +covMatrix[18]*2.*cp*ct*st // GetCovPxPz
3836  +covMatrix[14]*sp*sp*ct*ct // GetCovPyPy
3837  +covMatrix[19]*2.*sp*ct*st // GetCovPyPz
3838  +covMatrix[20]*st*st; // GetCovPzPz
3839  Double_t normDecayTimeToDStar = secondDaughter->NormalizedDecayLength(vertexDStar) / (299792458 * TMath::Sqrt(1/((pdgMassMother*pdgMassMother*errorMomentum*errorMomentum/(momentumMother*momentumMother)) + 1)));
3840 
3841  ((TH1F*)fMotherHistogramArray[motherType][histType][26])->Fill(pointingAngleToDStar);
3842  ((TH1F*)fMotherHistogramArray[motherType][histType][27])->Fill(d0D0);
3843  ((TH1F*)fMotherHistogramArray[motherType][histType][28])->Fill(d0FirstDaughterToDStar);
3844  ((TH1F*)fMotherHistogramArray[motherType][histType][29])->Fill(d0SecondDaughterToDStar);
3845  ((TH1F*)fMotherHistogramArray[motherType][histType][30])->Fill(impactProductToDStar);
3846  ((TH1F*)fMotherHistogramArray[motherType][histType][31])->Fill(impactProductXYToDStar);
3847  ((TH1F*)fMotherHistogramArray[motherType][histType][32])->Fill(normDecayLengthToDStar);
3848  ((TH1F*)fMotherHistogramArray[motherType][histType][33])->Fill(pseudoProperDecayTimeToDStar);
3849  ((TH1F*)fMotherHistogramArray[motherType][histType][34])->Fill(DecayTimeToDStar);
3850  ((TH1F*)fMotherHistogramArray[motherType][histType][35])->Fill(normDecayTimeToDStar);
3851 
3852 
3853  }
3854  return;
3855 }
3856 //-------------------------------------------------------------------------------------
3857 Int_t AliAnalysisTaskSEB0toDStarPi::MatchCandidateToMonteCarlo(Int_t pdgabs, AliAODRecoDecayHF2Prong * candidate, TClonesArray *mcArray, TMatrix * B0toDStarPiLabelMatrix) const
3858 {
3859  //
3860  // Check if this candidate is matched to a MC signal
3861  // If no, return -1
3862  // If yes, return label (>=0) of the AliAODMCParticle
3863 
3864 
3865  // Check number of daughters
3866  Int_t ndg = candidate->GetNDaughters();
3867  if(!ndg) { AliError("No daughters available"); return -1;}
3868  if(ndg != 2) return -1;
3869 
3870  // Skip if there is no signal in the event
3871  if(B0toDStarPiLabelMatrix->GetNrows() == 0) return -1;
3872 
3873  // loop on daughters and write the labels
3874  Int_t dgLabels[2] = {-1};
3875  Int_t pdgDg[2] = {0};
3876  Int_t signalPosition = -1;
3877  if(pdgabs==421)
3878  {
3879  AliAODTrack *trk0 = (AliAODTrack*)candidate->GetDaughter(0);
3880  dgLabels[0] = trk0->GetLabel();
3881  AliAODTrack *trk1 = (AliAODTrack*)candidate->GetDaughter(1);
3882  dgLabels[1] = trk1->GetLabel();
3883  pdgDg[0] = 211; pdgDg[1] = 321;
3884  signalPosition = 4;
3885  }
3886  else if(pdgabs==413)
3887  {
3888  AliAODTrack *trk0 = (AliAODTrack*)candidate->GetDaughter(0);
3889  dgLabels[0] = trk0->GetLabel();
3890  dgLabels[1] = MatchCandidateToMonteCarlo(421,(AliAODRecoDecayHF2Prong*)candidate->GetDaughter(1), mcArray, B0toDStarPiLabelMatrix);
3891  pdgDg[0] = 211; pdgDg[1] = 421;
3892  signalPosition = 5;
3893  }
3894  else if(pdgabs==511)
3895  {
3896  AliAODTrack *trk0 = (AliAODTrack*)candidate->GetDaughter(0);
3897  dgLabels[0] = trk0->GetLabel();
3898  dgLabels[1] = MatchCandidateToMonteCarlo(413,(AliAODRecoDecayHF2Prong*)candidate->GetDaughter(1), mcArray, B0toDStarPiLabelMatrix);
3899  pdgDg[0] = 211; pdgDg[1] = 413;
3900  signalPosition = 6;
3901  }
3902  else
3903  {
3904  std::cout << "Wrong pdg supplied for function to match candidate to monte carlo signal." << std::endl;
3905  return -1;
3906  }
3907  if(dgLabels[0]==-1) return -1;
3908  if(dgLabels[1]==-1) return -1;
3909 
3910 
3911  Int_t labMom[2]={0,0};
3912  Int_t i,j,lab,labMother,pdgMother,pdgPart;
3913  AliAODMCParticle *part=0;
3914  AliAODMCParticle *mother=0;
3915  Double_t pxSumDgs=0.,pySumDgs=0.,pzSumDgs=0.;
3916  Bool_t pdgUsed[2]={kFALSE,kFALSE};
3917 
3918  // loop on daughter labels
3919  for(i=0; i<ndg; i++)
3920  {
3921  labMom[i]=-1;
3922  lab = TMath::Abs(dgLabels[i]);
3923  if(lab<0)
3924  {
3925  printf("daughter with negative label %d\n",lab);
3926  return -1;
3927  }
3928  part = (AliAODMCParticle*)mcArray->At(lab);
3929  if(!part)
3930  {
3931  printf("no MC particle\n");
3932  return -1;
3933  }
3934 
3935  // check the PDG of the daughter
3936  pdgPart=TMath::Abs(part->GetPdgCode());
3937  for(j=0; j<ndg; j++)
3938  {
3939  if(!pdgUsed[j] && pdgPart==pdgDg[j])
3940  {
3941  pdgUsed[j]=kTRUE;
3942  break;
3943  }
3944  }
3945 
3946 
3947  mother = part;
3948  while(mother->GetMother()>=0)
3949  {
3950  labMother=mother->GetMother();
3951  mother = (AliAODMCParticle*)mcArray->At(labMother);
3952  if(!mother)
3953  {
3954  printf("no MC mother particle\n");
3955  break;
3956  }
3957  pdgMother = TMath::Abs(mother->GetPdgCode());
3958  if(pdgMother==pdgabs)
3959  {
3960  labMom[i]=labMother;
3961  // keep sum of daughters' momenta, to check for mom conservation
3962  pxSumDgs += part->Px();
3963  pySumDgs += part->Py();
3964  pzSumDgs += part->Pz();
3965  break;
3966  }
3967  else break;
3968  }
3969  if(labMom[i]==-1) return -1; // mother PDG not ok for this daughter
3970  } // end loop on daughters
3971 
3972  // check if the candidate is signal
3973  labMother=labMom[0];
3974  // all labels have to be the same and !=-1
3975  for(i=0; i<ndg; i++)
3976  {
3977  if(labMom[i]==-1) return -1;
3978  if(labMom[i]!=labMother) return -1;
3979  }
3980 
3981  // check that all daughter PDGs are matched
3982  for(i=0; i<ndg; i++)
3983  {
3984  if(pdgUsed[i]==kFALSE) return -1;
3985  }
3986 
3987  // Check for mom conservation
3988  mother = (AliAODMCParticle*)mcArray->At(labMother);
3989  Double_t pxMother = mother->Px();
3990  Double_t pyMother = mother->Py();
3991  Double_t pzMother = mother->Pz();
3992  // within 0.1%
3993  if((TMath::Abs(pxMother-pxSumDgs)/(TMath::Abs(pxMother)+1.e-13)) > 0.001 &&
3994  (TMath::Abs(pyMother-pySumDgs)/(TMath::Abs(pyMother)+1.e-13)) > 0.001 &&
3995  (TMath::Abs(pzMother-pzSumDgs)/(TMath::Abs(pzMother)+1.e-13)) > 0.001)
3996  {
3997  // cout << "Difference: " <<
3998  // (TMath::Abs(pxMother-pxSumDgs)/(TMath::Abs(pxMother)+1.e-13)) << " " <<
3999  // (TMath::Abs(pyMother-pySumDgs)/(TMath::Abs(pyMother)+1.e-13)) << " " <<
4000  // (TMath::Abs(pzMother-pzSumDgs)/(TMath::Abs(pzMother)+1.e-13)) << endl;
4001  return -1;
4002  }
4003 
4004  // Check if label matches a signal label
4005  Int_t bIsSignal = kFALSE;
4006  TMatrix &particleMatrix = *B0toDStarPiLabelMatrix;
4007  for (Int_t k = 0; k < B0toDStarPiLabelMatrix->GetNrows(); ++k)
4008  {
4009  if(labMother == (Int_t)particleMatrix(k,signalPosition))
4010  {
4011  bIsSignal = kTRUE;
4012  break;
4013  }
4014  }
4015  if(!bIsSignal) return -1;
4016 
4017  return labMother;
4018 }
Double_t NormalizedDecayLength() const
void FillD0Histograms(AliAODRecoDecayHF2Prong *selectedMother, AliAODVertex *primaryVertex, Double_t bz, Int_t motherType, Int_t histType, Int_t pdgCodeMother=-1)
virtual Int_t SelectPID(AliAODTrack *track, Int_t type)
double Double_t
Definition: External.C:58
ULong64_t GetTriggerMask()
Definition: AliRDHFCuts.h:69
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 FillDStarAndB0Histograms(AliAODRecoDecayHF2Prong *selectedMother, AliAODVertex *primaryVertex, Double_t bz, Int_t motherType, Int_t histType)
virtual void Terminate(Option_t *option)
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
Int_t ApplyCutOnVariableDStarforDStarptbin(Int_t nCutIndex, Int_t ptbin, Float_t cutVariableValue, Bool_t bCutArray[25])
Int_t MatchCandidateToMonteCarlo(Int_t pdgabs, AliAODRecoDecayHF2Prong *candidate, TClonesArray *mcArray, TMatrix *B0toDStarPiLabelMatrix) const
void InvMass(Int_t icalo, TString particle, TString fileName)
AliNormalizationCounter * fCounter
!Counter for normalization slot 4
Double_t DeltaInvMassB0Kpipipi(AliAODRecoDecayHF2Prong *B0) const
Int_t PtBinDStarforDStarptbin(Double_t pt) const
Int_t GetWhyRejection() const
Definition: AliRDHFCuts.h:313
Int_t IsD0forD0ptbinSelected(TObject *obj, Int_t selectionLevel, AliAODEvent *aod, Bool_t *bCutArray)
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
virtual void UserCreateOutputObjects()
Implementation of interface methods.
Int_t PtBinD0forDStarptbin(Double_t pt) const
TString part
use mixed event to constrain combinatorial background
Definition: InvMassFit.C:52
AliAODPidHF * GetPidHF() const
Definition: AliRDHFCuts.h:246
void D0Selection(AliAODEvent *aodEvent, AliAODVertex *primaryVertex, Double_t bz, TClonesArray *mcTrackArray, TMatrix *B0toDStarPiLabelMatrix, TClonesArray *D0TracksFromFriendFile)
int Int_t
Definition: External.C:63
Definition: External.C:204
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel, AliAODEvent *aod, Bool_t *bCutArray)
void B0toDStarPiSignalTracksInMC(TClonesArray *mcTrackArray, AliAODEvent *aodevent, TMatrix *B0toDStarPiLabelMatrix, TList *listout)
Int_t IsDStarforDStarptbinSelected(TObject *obj, Int_t selectionLevel, AliAODEvent *aod, Bool_t *bCutArray)
Float_t * GetPtBinLimitsDStarforDStarptbin() const
Int_t GetNPtBinsD0forD0ptbin() const
void DStarAndB0Selection(AliAODEvent *aodEvent, AliAODVertex *primaryVertex, Double_t bz, TClonesArray *mcTrackArray, TMatrix *B0toDStarPiLabelMatrix, TClonesArray *D0TracksFromFriendFile)
Float_t * GetPtBinLimitsD0forD0ptbin() const
void SetProngIDs(Int_t nIDs, UShort_t *id)
UShort_t GetProngID(Int_t ip) const
virtual void UserExec(Option_t *option)
void SetPrimaryVtxRef(TObject *vtx)
primary vertex
short Short_t
Definition: External.C:23
void B0PionSelection(AliAODEvent *aodEvent, AliAODVertex *primaryVertex, Double_t bz, TClonesArray *mcTrackArray, TMatrix *B0toDStarPiLabelMatrix)
void FillFinalTrackHistograms(AliAODRecoDecayHF2Prong *trackB0, Bool_t isDesiredCandidate, TClonesArray *mcTrackArray)
Bool_t IsEventSelected(AliVEvent *event)
AliAODVertex * RecalculateVertex(const AliVVertex *primary, TObjArray *tracks, Double_t bField, Double_t dispersion)
void StoreEvent(AliVEvent *, AliRDHFCuts *, Bool_t mc=kFALSE, Int_t multiplicity=-9999, Double_t spherocity=-99.)
Int_t PtBinD0forD0ptbin(Double_t pt) const
TH1F * fCEvents
Cuts - sent to output slot 3.
Float_t * GetPtBinLimits() const
Definition: AliRDHFCuts.h:247
Int_t GetNPtBinsDStarforDStarptbin() const
Double_t DeltaInvMassDStarKpipi(AliAODRecoDecayHF2Prong *DStar) const
unsigned short UShort_t
Definition: External.C:28
Int_t IsD0forDStarptbinSelected(TObject *obj, Int_t selectionLevel, AliAODEvent *aod, Bool_t *bCutArray)
const char Option_t
Definition: External.C:48
Float_t * GetPtBinLimitsD0forDStarptbin() const
Int_t GetNPtBinsD0forDStarptbin() const
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:248
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
Int_t PtBin(Double_t pt) const
void DStarPionSelection(AliAODEvent *aodEvent, AliAODVertex *primaryVertex, Double_t bz, TClonesArray *mcTrackArray, TMatrix *B0toDStarPiLabelMatrix)