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