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