AliPhysics  a6017e1 (a6017e1)
AliAnalysisTaskSEDStarSpectra.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 DStar Analysis
21 //
22 //
23 // The D* spectra study is done in pt bins:
24 // [0,0.5] [0.5,1] [1,2] [2,3] [3,4] [4,5] [5,6] [6,7] [7,8],
25 // [8,10],[10,12], [12,16], [16,20] and [20,24]
26 //
27 // Cuts arew centralized in AliRDHFCutsDStartoKpipi
28 // Side Band and like sign background are implemented in the macro
29 //
30 //-----------------------------------------------------------------------
31 //
32 // Author A.Grelli
33 // ERC-QGP Utrecht University - a.grelli@uu.nl,
34 // Author Y.Wang
35 // University of Heidelberg - yifei@physi.uni-heidelberg.de
36 // Author C.Ivan
37 // ERC-QGP Utrecht University - c.ivan@uu.nl,
38 //
39 //-----------------------------------------------------------------------
40 
41 #include <TSystem.h>
42 #include <TParticle.h>
43 #include <TH1I.h>
44 #include "TROOT.h"
45 #include <TDatabasePDG.h>
46 #include <AliAnalysisDataSlot.h>
47 #include <AliAnalysisDataContainer.h>
49 #include "AliMCEvent.h"
50 #include "AliAnalysisManager.h"
51 #include "AliAODMCHeader.h"
52 #include "AliAODHandler.h"
53 #include "AliLog.h"
54 #include "AliAODVertex.h"
55 #include "AliAODRecoDecay.h"
56 #include "AliAODRecoDecayHF.h"
57 #include "AliAODRecoCascadeHF.h"
59 #include "AliAnalysisVertexingHF.h"
60 #include "AliESDtrack.h"
61 #include "AliAODMCParticle.h"
63 #include "AliAODEvent.h"
65 
69 
70 //__________________________________________________________________________
73  fEvents(0),
74  fAnalysis(0),
75  fD0Window(0),
76  fPeakWindow(0),
77  fUseMCInfo(kFALSE),
78  fDoSearch(kFALSE),
79  fOutput(0),
80  fOutputAll(0),
81  fOutputPID(0),
82  fNSigma(3),
83  fCuts(0),
84  fCEvents(0),
85  fTrueDiff2(0),
86  fDeltaMassD1(0),
87  fCounter(0),
88  fAODProtection(1),
89  fDoImpParDstar(kFALSE),
90  fNImpParBins(400),
91  fLowerImpPar(-2000.),
92  fHigherImpPar(2000.),
93  fNPtBins(0),
94  fAllhist(0x0),
95  fPIDhist(0x0),
96  fDoDStarVsY(kFALSE)
97 {
98  //
100  //
101  for(Int_t i=0;i<5;i++) fHistMassPtImpParTCDs[i]=0;
102 }
103 
104 //___________________________________________________________________________
106  AliAnalysisTaskSE(name),
107  fEvents(0),
108  fAnalysis(0),
109  fD0Window(0),
110  fPeakWindow(0),
111  fUseMCInfo(kFALSE),
112  fDoSearch(kFALSE),
113  fOutput(0),
114  fOutputAll(0),
115  fOutputPID(0),
116  fNSigma(3),
117  fCuts(0),
118  fCEvents(0),
119  fTrueDiff2(0),
120  fDeltaMassD1(0),
121  fCounter(0),
122  fAODProtection(1),
123  fDoImpParDstar(kFALSE),
124  fNImpParBins(400),
125  fLowerImpPar(-2000.),
126  fHigherImpPar(2000.),
127  fNPtBins(0),
128  fAllhist(0x0),
129  fPIDhist(0x0),
130  fDoDStarVsY(kFALSE)
131 {
132  //
134  //
135  Info("AliAnalysisTaskSEDStarSpectra","Calling Constructor");
136 
137  fCuts=cuts;
138  for(Int_t i=0;i<5;i++) fHistMassPtImpParTCDs[i]=0;
139 
140 
141  DefineOutput(1,TList::Class()); //counters
142  DefineOutput(2,TList::Class()); //All Entries output
143  DefineOutput(3,TList::Class()); //3sigma PID output
144  DefineOutput(4,AliRDHFCutsDStartoKpipi::Class()); //My private output
145  DefineOutput(5,AliNormalizationCounter::Class()); // normalization
146 }
147 
148 //___________________________________________________________________________
150  //
152  //
153  Info("~AliAnalysisTaskSEDStarSpectra","Calling Destructor");
154 
155  delete fOutput;
156  delete fOutputAll;
157  delete fOutputPID;
158  delete fCuts;
159  delete fCEvents;
160  delete fDeltaMassD1;
161  for(Int_t i=0; i<5; i++){
162  delete fHistMassPtImpParTCDs[i];
163  }
164  for(Int_t i=0; i<((fNPtBins+2)*18); i++){
165  delete fAllhist[i];
166  delete fPIDhist[i];
167  }
168  delete [] fAllhist;
169  delete [] fPIDhist;
170 
171 }
172 //_________________________________________________
174  //
176  //
177 
178  if(fDebug > 1) printf("AnalysisTaskSEDStarSpectra::Init() \n");
181  // Post the data
182  PostData(4,copyfCuts);
183 
184  return;
185 }
186 
187 //_________________________________________________
189 {
191  if (!fInputEvent) {
192  Error("UserExec","NO EVENT FOUND!");
193  return;
194  }
195 
196 
197  if(fAODProtection>=0){
198  // Protection against different number of events in the AOD and deltaAOD
199  // In case of discrepancy the event is rejected.
200  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
201  if (matchingAODdeltaAODlevel<0 || (matchingAODdeltaAODlevel==0 && fAODProtection==1)) {
202  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
203  fCEvents->Fill(11);
204  return;
205  }
206  }
207 
208  fEvents++;
209 
210  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
211  TClonesArray *arrayDStartoD0pi=0;
212 
213  fCEvents->Fill(1);
214 
215  if(!aodEvent && AODEvent() && IsStandardAOD()) {
216  // In case there is an AOD handler writing a standard AOD, use the AOD
217  // event in memory rather than the input (ESD) event.
218  aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
219  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
220  // have to taken from the AOD event hold by the AliAODExtension
221  AliAODHandler* aodHandler = (AliAODHandler*)
222  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
223  if(aodHandler->GetExtensions()) {
224  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
225  AliAODEvent *aodFromExt = ext->GetAOD();
226  arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
227  }
228  } else {
229  arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
230  }
231 
232  // fix for temporary bug in ESDfilter
233  // the AODs with null vertex pointer didn't pass the PhysSel
234  if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
235  fCEvents->Fill(2);
236 
237  fCounter->StoreEvent(aodEvent,fCuts,fUseMCInfo);
238 
239  // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
240  TString trigclass=aodEvent->GetFiredTriggerClasses();
241  if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")||trigclass.Contains("C0SMH-B-NOPF-ALL")) fCEvents->Fill(5);
242 
243  if(!fCuts->IsEventSelected(aodEvent)) {
244  if(fCuts->GetWhyRejection()==6) // rejected for Z vertex
245  fCEvents->Fill(6);
246  return;
247  }
248 
249  Bool_t isEvSel=fCuts->IsEventSelected(aodEvent);
250  fCEvents->Fill(3);
251  if(!isEvSel) return;
252 
253  // Load the event
254  // AliInfo(Form("Event %d",fEvents));
255  //if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
256 
257  // counters for efficiencies
258  Int_t icountReco = 0;
259 
260  //D* and D0 prongs needed to MatchToMC method
261  Int_t pdgDgDStartoD0pi[2]={421,211};
262  Int_t pdgDgD0toKpi[2]={321,211};
263 
264  // AOD primary vertex
265  AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
266  if(!vtx1) return;
267  if(vtx1->GetNContributors()<1) return;
268  fCEvents->Fill(4);
269 
270  if (!arrayDStartoD0pi){
271  AliInfo("Could not find array of HF vertices, skipping the event");
272  return;
273  }else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
274 
275  Int_t nSelectedAna =0;
276  Int_t nSelectedProd =0;
277 
278  // vHF object is needed to call the method that refills the missing info of the candidates
279  // if they have been deleted in dAOD reconstruction phase
280  // in order to reduce the size of the file
282 
283  // loop over the tracks to search for candidates soft pion
284  for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
285 
286  // D* candidates and D0 from D*
287  AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
288  Bool_t isDStarCand =kTRUE;
289  if(!(vHF->FillRecoCasc(aodEvent,dstarD0pi,isDStarCand))) {//Fill the data members of the candidate only if they are empty.
290  fCEvents->Fill(12); //monitor how often this fails
291  continue;
292  }
293  if(!dstarD0pi->GetSecondaryVtx()) continue;
294  AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
295  if (!theD0particle) continue;
296 
297  Int_t isDStar = 0;
298  TClonesArray *mcArray = 0;
299  AliAODMCHeader *mcHeader=0;
300 
301  Bool_t isPrimary=kTRUE;
302  Float_t pdgCode=-2;
303  Float_t trueImpParXY=0.;
304 
305  // mc analysis
306  if(fUseMCInfo){
307  //MC array need for maching
308  mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
309  if (!mcArray) {
310  AliError("Could not find Monte-Carlo in AOD");
311  return;
312  }
313  // load MC header
314  mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
315  if(!mcHeader) {
316  printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
317  return;
318  }
319  // find associated MC particle for D* ->D0toKpi
320  Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
321  if(mcLabel>=0){
322 
323  AliAODMCParticle *partDSt = (AliAODMCParticle*)mcArray->At(mcLabel);
324  Int_t checkOrigin = CheckOrigin(mcArray,partDSt);
325  if(checkOrigin==5) isPrimary=kFALSE;
326  AliAODMCParticle *dg0 = (AliAODMCParticle*)mcArray->At(partDSt->GetDaughter(0));
327  // AliAODMCParticle *dg01 = (AliAODMCParticle*)mcArray->At(dg0->GetDaughter(0));
328  pdgCode=TMath::Abs(partDSt->GetPdgCode());
329  if(!isPrimary){
330  trueImpParXY=GetTrueImpactParameterD0(mcHeader,mcArray,dg0)*1000.;
331  }
332  isDStar = 1;
333  }else{
334  pdgCode=-1;
335  }
336  }
337 
338  if(pdgCode==-1) AliDebug(2,"No particle assigned! check\n");
339 
340  Double_t Dstarpt = dstarD0pi->Pt();
341 
342  // quality selction on tracks and region of interest
343  Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
344  if(!isTkSelected) continue;
345 
346  if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
347 
348 
349  //histos for impact par studies - D0!!!
350  Double_t ptCand = dstarD0pi->Get2Prong()->Pt();
351  Double_t invMass=dstarD0pi->InvMassD0();
352  Double_t impparXY=dstarD0pi->Get2Prong()->ImpParXY()*10000.;
353 
354  Double_t arrayForSparse[3]={invMass,ptCand,impparXY};
355  Double_t arrayForSparseTrue[3]={invMass,ptCand,trueImpParXY};
356 
357  // set the D0 and D* search window bin by bin - D* window useful to speed up the reconstruction and D0 window used *ONLY* to calculate side band bkg for the background subtraction methods, for the standard analysis the value in the cut file is considered
358 
359  if (0<=Dstarpt && Dstarpt<0.5){
360  if(fAnalysis==1){
361  fD0Window=0.035;
362  fPeakWindow=0.03;
363  }else{
364  fD0Window=0.020;
365  fPeakWindow=0.0018;
366  }
367  }
368  if (0.5<=Dstarpt && Dstarpt<1.0){
369  if(fAnalysis==1){
370  fD0Window=0.035;
371  fPeakWindow=0.03;
372  }else{
373  fD0Window=0.020;
374  fPeakWindow=0.0018;
375  }
376  }
377  if (1.0<=Dstarpt && Dstarpt<2.0){
378  if(fAnalysis==1){
379  fD0Window=0.035;
380  fPeakWindow=0.03;
381  }else{
382  fD0Window=0.020;
383  fPeakWindow=0.0018;
384  }
385  }
386  if (2.0<=Dstarpt && Dstarpt<3.0){
387  if(fAnalysis==1){
388  fD0Window=0.035;
389  fPeakWindow=0.03;
390  }else{
391  fD0Window=0.022;
392  fPeakWindow=0.0016;
393  }
394  }
395  if (3.0<=Dstarpt && Dstarpt<4.0){
396  if(fAnalysis==1){
397  fD0Window=0.035;
398  fPeakWindow=0.03;
399  }else{
400  fD0Window=0.026;
401  fPeakWindow=0.0014;
402  }
403  }
404  if (4.0<=Dstarpt && Dstarpt<5.0){
405  if(fAnalysis==1){
406  fD0Window=0.045;
407  fPeakWindow=0.03;
408  }else{
409  fD0Window=0.026;
410  fPeakWindow=0.0014;
411  }
412  }
413  if (5.0<=Dstarpt && Dstarpt<6.0){
414  if(fAnalysis==1){
415  fD0Window=0.045;
416  fPeakWindow=0.03;
417  }else{
418  fD0Window=0.026;
419  fPeakWindow=0.006;
420  }
421  }
422  if (6.0<=Dstarpt && Dstarpt<7.0){
423  if(fAnalysis==1){
424  fD0Window=0.055;
425  fPeakWindow=0.03;
426  }else{
427  fD0Window=0.026;
428  fPeakWindow=0.006;
429  }
430  }
431  if (Dstarpt>=7.0){
432  if(fAnalysis==1){
433  fD0Window=0.074;
434  fPeakWindow=0.03;
435  }else{
436  fD0Window=0.026;
437  fPeakWindow=0.006;
438  }
439  }
440 
441  nSelectedProd++;
442  nSelectedAna++;
443 
444  // check that we are close to signal in the DeltaM - here to save time for PbPb
445  Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
446  Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
447  Double_t invmassDelta = dstarD0pi->DeltaInvMass();
448 
449  if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))>fPeakWindow) continue;
450  Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate,aodEvent); //selected
451 
452  // after cuts
453  if(fDoImpParDstar && isSelected){
454  fHistMassPtImpParTCDs[0]->Fill(arrayForSparse);
455  if(isPrimary) fHistMassPtImpParTCDs[1]->Fill(arrayForSparse);
456  else{
457  fHistMassPtImpParTCDs[2]->Fill(arrayForSparse);
458  fHistMassPtImpParTCDs[3]->Fill(arrayForSparseTrue);
459  }
460  }
461 
462  if (fDoDStarVsY && isSelected){
463  ((TH3F*) (fOutputPID->FindObject("deltamassVsyVsPt")))->Fill(dstarD0pi->DeltaInvMass(),dstarD0pi->YDstar(),dstarD0pi->Pt() );
464  }
465 
466 
467  // fill PID
468  FillSpectrum(dstarD0pi,isDStar,fCuts,isSelected,fOutputPID, fPIDhist);
469  SideBandBackground(dstarD0pi,fCuts,isSelected, fOutputPID, fPIDhist);
470  //WrongSignForDStar(dstarD0pi,fCuts,fOutputPID);
471 
472  //swich off the PID selection
473  fCuts->SetUsePID(kFALSE);
474  Int_t isSelectedNoPID=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate, aodEvent); //selected
475  fCuts->SetUsePID(kTRUE);
476 
477  FillSpectrum(dstarD0pi,isDStar,fCuts,isSelectedNoPID,fOutputAll, fAllhist);
478  // SideBandBackground(dstarD0pi,fCuts,isSelectedNoPID, fOutputAll);
479 
480  // rare D search ------
481  if(fDoSearch){
482  TLorentzVector lorentzTrack1(0,0,0,0); // lorentz 4 vector
483  TLorentzVector lorentzTrack2(0,0,0,0); // lorentz 4 vector
484 
485  for (Int_t i=0; i<aodEvent->GetNumberOfTracks(); i++){
486 
487  AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(i));
488  if(!aodTrack) AliFatal("Not a standard AOD");
489 
490  if(dstarD0pi->Charge() == aodTrack->Charge()) continue;
491  if((!(aodTrack->GetStatus()&AliESDtrack::kITSrefit)|| (!(aodTrack->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
492  if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))>0.02) continue;
493 
494  //build the D1 mass
495  Double_t mass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
496 
497  lorentzTrack1.SetPxPyPzE( dstarD0pi->Px(),dstarD0pi->Py(), dstarD0pi->Pz(), dstarD0pi->E(413) );
498  lorentzTrack2.SetPxPyPzE( aodTrack->Px(),aodTrack->Py(), aodTrack->Pz(),aodTrack->E(mass) );
499 
500  //D1 mass
501  Double_t d1mass = ((lorentzTrack1+lorentzTrack2).M());
502  //mass difference - at 0.4117 and 0.4566
503  fDeltaMassD1->Fill(d1mass-dstarD0pi->InvMassDstarKpipi());
504  }
505  }
506 
507  if(isDStar == 1) {
508  fTrueDiff2->Fill(dstarD0pi->Pt(),dstarD0pi->DeltaInvMass());
509  }
510 
511  }
512 
513  fCounter->StoreCandidates(aodEvent,nSelectedProd,kTRUE);
514  fCounter->StoreCandidates(aodEvent,nSelectedAna,kFALSE);
515 
516  delete vHF;
517 
518  AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
519 
520  PostData(1,fOutput);
521  PostData(2,fOutputAll);
522  PostData(3,fOutputPID);
523  PostData(5,fCounter);
524 
525 }
526 //________________________________________ terminate ___________________________
528 {
532 
533  //Info("Terminate","");
534  AliAnalysisTaskSE::Terminate();
535 
536  fOutput = dynamic_cast<TList*> (GetOutputData(1));
537  if (!fOutput) {
538  printf("ERROR: fOutput not available\n");
539  return;
540  }
541 
542  fCEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fCEvents"));
543  fDeltaMassD1 = dynamic_cast<TH1F*>(fOutput->FindObject("fDeltaMassD1"));
544  fTrueDiff2 = dynamic_cast<TH2F*>(fOutput->FindObject("fTrueDiff2"));
545 
546  fOutputAll = dynamic_cast<TList*> (GetOutputData(1));
547  if (!fOutputAll) {
548  printf("ERROR: fOutputAll not available\n");
549  return;
550  }
551  fOutputPID = dynamic_cast<TList*> (GetOutputData(2));
552  if (!fOutputPID) {
553  printf("ERROR: fOutputPID not available\n");
554  return;
555  }
556  return;
557 }
558 //___________________________________________________________________________
561  Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
562 
563  //slot #1
564  //OpenFile(1);
565  fOutput = new TList();
566  fOutput->SetOwner();
567  fOutput->SetName("chist0");
568 
569  fOutputAll = new TList();
570  fOutputAll->SetOwner();
571  fOutputAll->SetName("listAll");
572 
573  fOutputPID = new TList();
574  fOutputPID->SetOwner();
575  fOutputPID->SetName("listPID");
576 
577  // define histograms
579 
580  //Counter for Normalization
581  fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));
582  fCounter->Init();
583 
585 
586  PostData(1,fOutput);
587  PostData(2,fOutputAll);
588  PostData(3,fOutputPID);
589 
590  return;
591 }
592 //___________________________________ hiostograms _______________________________________
595 
596  fCEvents = new TH1F("fCEvents","counter",13,0,13);
597  fCEvents->SetStats(kTRUE);
598  fCEvents->GetXaxis()->SetTitle("1");
599  fCEvents->GetYaxis()->SetTitle("counts");
600  fCEvents->GetXaxis()->SetBinLabel(2,"no. of events");
601  fCEvents->GetXaxis()->SetBinLabel(3,"good prim vtx and B field");
602  fCEvents->GetXaxis()->SetBinLabel(4,"no event selected");
603  fCEvents->GetXaxis()->SetBinLabel(5,"no vtx contributors");
604  fCEvents->GetXaxis()->SetBinLabel(6,"trigger for PbPb");
605  fCEvents->GetXaxis()->SetBinLabel(7,"no z vtx");
606  fCEvents->GetXaxis()->SetBinLabel(12,"no. of D0 fail to be rec");
607  fOutput->Add(fCEvents);
608 
609  fTrueDiff2 = new TH2F("DiffDstar_pt","True Reco diff vs pt",200,0,15,900,0,0.3);
610  fOutput->Add(fTrueDiff2);
611 
612  fDeltaMassD1 = new TH1F("DeltaMassD1","delta mass d1",600,0,0.8);
613  fOutput->Add(fDeltaMassD1);
614  //temp a
615 
616  fAllhist = new TH1F*[(fNPtBins+2)*18];
617  fPIDhist = new TH1F*[(fNPtBins+2)*18];
618 
619  TString nameMass=" ", nameSgn=" ", nameBkg=" ";
620 
621  for(Int_t i=-2;i<fNPtBins;i++){
622  nameMass="histDeltaMass_";
623  nameMass+=i+1;
624  nameSgn="histDeltaSgn_";
625  nameSgn+=i+1;
626  nameBkg="histDeltaBkg_";
627  nameBkg+=i+1;
628 
629  if (i==-2) {
630  nameMass="histDeltaMass";
631  nameSgn="histDeltaSgn";
632  nameBkg="histDeltaBkg";
633  }
634 
635  TH1F* spectrumMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} invariant mass; #DeltaM [GeV/c^{2}]; Entries",700,0.13,0.2);
636  TH1F* spectrumSgn = new TH1F(nameSgn.Data(), "D^{*}-D^{0} Signal invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",700,0.13,0.2);
637  TH1F* spectrumBkg = new TH1F(nameBkg.Data(), "D^{*}-D^{0} Background invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",700,0.13,0.2);
638 
639  nameMass="histD0Mass_";
640  nameMass+=i+1;
641  nameSgn="histD0Sgn_";
642  nameSgn+=i+1;
643  nameBkg="histD0Bkg_";
644  nameBkg+=i+1;
645 
646  if (i==-2) {
647  nameMass="histD0Mass";
648  nameSgn="histD0Sgn";
649  nameBkg="histD0Bkg";
650  }
651 
652  TH1F* spectrumD0Mass = new TH1F(nameMass.Data(),"D^{0} invariant mass; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
653  TH1F* spectrumD0Sgn = new TH1F(nameSgn.Data(), "D^{0} Signal invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
654  TH1F* spectrumD0Bkg = new TH1F(nameBkg.Data(), "D^{0} Background invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
655 
656  nameMass="histDstarMass_";
657  nameMass+=i+1;
658  nameSgn="histDstarSgn_";
659  nameSgn+=i+1;
660  nameBkg="histDstarBkg_";
661  nameBkg+=i+1;
662 
663  if (i==-2) {
664  nameMass="histDstarMass";
665  nameSgn="histDstarSgn";
666  nameBkg="histDstarBkg";
667  }
668 
669  TH1F* spectrumDstarMass = new TH1F(nameMass.Data(),"D^{*} invariant mass; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
670  TH1F* spectrumDstarSgn = new TH1F(nameSgn.Data(), "D^{*} Signal invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
671  TH1F* spectrumDstarBkg = new TH1F(nameBkg.Data(), "D^{*} Background invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
672 
673  nameMass="histSideBandMass_";
674  nameMass+=i+1;
675  if (i==-2) {
676  nameMass="histSideBandMass";
677  }
678 
679  TH1F* spectrumSideBandMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} sideband mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
680 
681  nameMass="histWrongSignMass_";
682  nameMass+=i+1;
683  if (i==-2) {
684  nameMass="histWrongSignMass";
685  }
686 
687  TH1F* spectrumWrongSignMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} wrongsign mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
688 
689 
690  spectrumMass->Sumw2();
691  spectrumSgn->Sumw2();
692  spectrumBkg->Sumw2();
693 
694  spectrumMass->SetLineColor(6);
695  spectrumSgn->SetLineColor(2);
696  spectrumBkg->SetLineColor(4);
697 
698  spectrumMass->SetMarkerStyle(20);
699  spectrumSgn->SetMarkerStyle(20);
700  spectrumBkg->SetMarkerStyle(20);
701  spectrumMass->SetMarkerSize(0.6);
702  spectrumSgn->SetMarkerSize(0.6);
703  spectrumBkg->SetMarkerSize(0.6);
704  spectrumMass->SetMarkerColor(6);
705  spectrumSgn->SetMarkerColor(2);
706  spectrumBkg->SetMarkerColor(4);
707 
708  spectrumD0Mass->Sumw2();
709  spectrumD0Sgn->Sumw2();
710  spectrumD0Bkg->Sumw2();
711 
712  spectrumD0Mass->SetLineColor(6);
713  spectrumD0Sgn->SetLineColor(2);
714  spectrumD0Bkg->SetLineColor(4);
715 
716  spectrumD0Mass->SetMarkerStyle(20);
717  spectrumD0Sgn->SetMarkerStyle(20);
718  spectrumD0Bkg->SetMarkerStyle(20);
719  spectrumD0Mass->SetMarkerSize(0.6);
720  spectrumD0Sgn->SetMarkerSize(0.6);
721  spectrumD0Bkg->SetMarkerSize(0.6);
722  spectrumD0Mass->SetMarkerColor(6);
723  spectrumD0Sgn->SetMarkerColor(2);
724  spectrumD0Bkg->SetMarkerColor(4);
725 
726  spectrumDstarMass->Sumw2();
727  spectrumDstarSgn->Sumw2();
728  spectrumDstarBkg->Sumw2();
729 
730  spectrumDstarMass->SetLineColor(6);
731  spectrumDstarSgn->SetLineColor(2);
732  spectrumDstarBkg->SetLineColor(4);
733 
734  spectrumDstarMass->SetMarkerStyle(20);
735  spectrumDstarSgn->SetMarkerStyle(20);
736  spectrumDstarBkg->SetMarkerStyle(20);
737  spectrumDstarMass->SetMarkerSize(0.6);
738  spectrumDstarSgn->SetMarkerSize(0.6);
739  spectrumDstarBkg->SetMarkerSize(0.6);
740  spectrumDstarMass->SetMarkerColor(6);
741  spectrumDstarSgn->SetMarkerColor(2);
742  spectrumDstarBkg->SetMarkerColor(4);
743 
744  spectrumSideBandMass->Sumw2();
745  spectrumSideBandMass->SetLineColor(4);
746  spectrumSideBandMass->SetMarkerStyle(20);
747  spectrumSideBandMass->SetMarkerSize(0.6);
748  spectrumSideBandMass->SetMarkerColor(4);
749 
750  spectrumWrongSignMass->Sumw2();
751  spectrumWrongSignMass->SetLineColor(4);
752  spectrumWrongSignMass->SetMarkerStyle(20);
753  spectrumWrongSignMass->SetMarkerSize(0.6);
754  spectrumWrongSignMass->SetMarkerColor(4);
755 
756  TH1F* allMass = (TH1F*)spectrumMass->Clone();
757  TH1F* allSgn = (TH1F*)spectrumSgn->Clone();
758  TH1F* allBkg = (TH1F*)spectrumBkg->Clone();
759 
760  TH1F* pidMass = (TH1F*)spectrumMass->Clone();
761  TH1F* pidSgn = (TH1F*)spectrumSgn->Clone();
762  TH1F* pidBkg = (TH1F*)spectrumBkg->Clone();
763 
764  fOutputAll->Add(allMass);
765  fOutputAll->Add(allSgn);
766  fOutputAll->Add(allBkg);
767  fAllhist[i+2+((fNPtBins+2)*kDeltaMass)]=allMass;
768  fAllhist[i+2+((fNPtBins+2)*kDeltaSgn)]=allSgn;
769  fAllhist[i+2+((fNPtBins+2)*kDeltaBkg)]=allBkg;
770 
771  fOutputPID->Add(pidMass);
772  fOutputPID->Add(pidSgn);
773  fOutputPID->Add(pidBkg);
774  fPIDhist[i+2+((fNPtBins+2)*kDeltaMass)]=pidMass;
775  fPIDhist[i+2+((fNPtBins+2)*kDeltaSgn)]=pidSgn;
776  fPIDhist[i+2+((fNPtBins+2)*kDeltaBkg)]=pidBkg;
777 
778  TH1F* allD0Mass = (TH1F*)spectrumD0Mass->Clone();
779  TH1F* allD0Sgn = (TH1F*)spectrumD0Sgn->Clone();
780  TH1F* allD0Bkg = (TH1F*)spectrumD0Bkg->Clone();
781 
782  TH1F* pidD0Mass = (TH1F*)spectrumD0Mass->Clone();
783  TH1F* pidD0Sgn = (TH1F*)spectrumD0Sgn->Clone();
784  TH1F* pidD0Bkg = (TH1F*)spectrumD0Bkg->Clone();
785 
786  fOutputAll->Add(allD0Mass);
787  fOutputAll->Add(allD0Sgn);
788  fOutputAll->Add(allD0Bkg);
789  fAllhist[i+2+((fNPtBins+2)*kDzMass)]=allD0Mass;
790  fAllhist[i+2+((fNPtBins+2)*kDzSgn)]=allD0Sgn;
791  fAllhist[i+2+((fNPtBins+2)*kDzBkg)]=allD0Bkg;
792 
793  fOutputPID->Add(pidD0Mass);
794  fOutputPID->Add(pidD0Sgn);
795  fOutputPID->Add(pidD0Bkg);
796  fPIDhist[i+2+((fNPtBins+2)*kDzMass)]=pidD0Mass;
797  fPIDhist[i+2+((fNPtBins+2)*kDzSgn)]=pidD0Sgn;
798  fPIDhist[i+2+((fNPtBins+2)*kDzBkg)]=pidD0Bkg;
799 
800  TH1F* allDstarMass = (TH1F*)spectrumDstarMass->Clone();
801  TH1F* allDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
802  TH1F* allDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
803 
804  TH1F* pidDstarMass = (TH1F*)spectrumDstarMass->Clone();
805  TH1F* pidDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
806  TH1F* pidDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
807 
808  fOutputAll->Add(allDstarMass);
809  fOutputAll->Add(allDstarSgn);
810  fOutputAll->Add(allDstarBkg);
811  fAllhist[i+2+((fNPtBins+2)*kDstarMass)]=allDstarMass;
812  fAllhist[i+2+((fNPtBins+2)*kDstarSgn)]=allDstarSgn;
813  fAllhist[i+2+((fNPtBins+2)*kDstarBkg)]=allDstarBkg;
814 
815  fOutputPID->Add(pidDstarMass);
816  fOutputPID->Add(pidDstarSgn);
817  fOutputPID->Add(pidDstarBkg);
818  fPIDhist[i+2+((fNPtBins+2)*kDstarMass)]=pidDstarMass;
819  fPIDhist[i+2+((fNPtBins+2)*kDstarSgn)]=pidDstarSgn;
820  fPIDhist[i+2+((fNPtBins+2)*kDstarBkg)]=pidDstarBkg;
821 
822  TH1F* allSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
823  TH1F* pidSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
824 
825  fOutputAll->Add(allSideBandMass);
826  fOutputPID->Add(pidSideBandMass);
827  fAllhist[i+2+((fNPtBins+2)*kSideBandMass)]=allSideBandMass;
828  fPIDhist[i+2+((fNPtBins+2)*kSideBandMass)]=pidSideBandMass;
829 
830  TH1F* allWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
831  TH1F* pidWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
832 
833  fOutputAll->Add(allWrongSignMass);
834  fOutputPID->Add(pidWrongSignMass);
835  fAllhist[i+2+((fNPtBins+2)*kWrongSignMass)]=allWrongSignMass;
836  fPIDhist[i+2+((fNPtBins+2)*kWrongSignMass)]=pidWrongSignMass;
837 
838  }
839 
840  // pt spectra
841  nameMass="ptMass";
842  nameSgn="ptSgn";
843  nameBkg="ptBkg";
844 
845  TH1F* ptspectrumMass = new TH1F(nameMass.Data(),"D^{*} p_{T}; p_{T} [GeV]; Entries",400,0,50);
846  TH1F* ptspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal p_{T} - MC; p_{T} [GeV]; Entries",400,0,50);
847  TH1F* ptspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background p_{T} - MC; p_{T} [GeV]; Entries",400,0,50);
848 
849  ptspectrumMass->Sumw2();
850  ptspectrumSgn->Sumw2();
851  ptspectrumBkg->Sumw2();
852 
853  ptspectrumMass->SetLineColor(6);
854  ptspectrumSgn->SetLineColor(2);
855  ptspectrumBkg->SetLineColor(4);
856 
857  ptspectrumMass->SetMarkerStyle(20);
858  ptspectrumSgn->SetMarkerStyle(20);
859  ptspectrumBkg->SetMarkerStyle(20);
860  ptspectrumMass->SetMarkerSize(0.6);
861  ptspectrumSgn->SetMarkerSize(0.6);
862  ptspectrumBkg->SetMarkerSize(0.6);
863  ptspectrumMass->SetMarkerColor(6);
864  ptspectrumSgn->SetMarkerColor(2);
865  ptspectrumBkg->SetMarkerColor(4);
866 
867  TH1F* ptallMass = (TH1F*)ptspectrumMass->Clone();
868  TH1F* ptallSgn = (TH1F*)ptspectrumSgn->Clone();
869  TH1F* ptallBkg = (TH1F*)ptspectrumBkg->Clone();
870 
871  TH1F* ptpidMass = (TH1F*)ptspectrumMass->Clone();
872  TH1F* ptpidSgn = (TH1F*)ptspectrumSgn->Clone();
873  TH1F* ptpidBkg = (TH1F*)ptspectrumBkg->Clone();
874 
875  fOutputAll->Add(ptallMass);
876  fOutputAll->Add(ptallSgn);
877  fOutputAll->Add(ptallBkg);
878  fAllhist[((fNPtBins+2)*kptMass)]=ptallMass;
879  fAllhist[((fNPtBins+2)*kptSgn)]=ptallSgn;
880  fAllhist[((fNPtBins+2)*kptBkg)]=ptallBkg;
881 
882  fOutputPID->Add(ptpidMass);
883  fOutputPID->Add(ptpidSgn);
884  fOutputPID->Add(ptpidBkg);
885  fPIDhist[(fNPtBins+2)*kptMass]=ptpidMass;
886  fPIDhist[(fNPtBins+2)*kptSgn]=ptpidSgn;
887  fPIDhist[(fNPtBins+2)*kptBkg]=ptpidBkg;
888  // eta spectra
889  nameMass="etaMass";
890  nameSgn="etaSgn";
891  nameBkg="etaBkg";
892 
893  TH1F* etaspectrumMass = new TH1F(nameMass.Data(),"D^{*} #eta; #eta; Entries",200,-1,1);
894  TH1F* etaspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal #eta - MC; #eta; Entries",200,-1,1);
895  TH1F* etaspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background #eta - MC; #eta; Entries",200,-1,1);
896 
897  etaspectrumMass->Sumw2();
898  etaspectrumSgn->Sumw2();
899  etaspectrumBkg->Sumw2();
900 
901  etaspectrumMass->SetLineColor(6);
902  etaspectrumSgn->SetLineColor(2);
903  etaspectrumBkg->SetLineColor(4);
904 
905  etaspectrumMass->SetMarkerStyle(20);
906  etaspectrumSgn->SetMarkerStyle(20);
907  etaspectrumBkg->SetMarkerStyle(20);
908  etaspectrumMass->SetMarkerSize(0.6);
909  etaspectrumSgn->SetMarkerSize(0.6);
910  etaspectrumBkg->SetMarkerSize(0.6);
911  etaspectrumMass->SetMarkerColor(6);
912  etaspectrumSgn->SetMarkerColor(2);
913  etaspectrumBkg->SetMarkerColor(4);
914 
915  TH1F* etaallMass = (TH1F*)etaspectrumMass->Clone();
916  TH1F* etaallSgn = (TH1F*)etaspectrumSgn->Clone();
917  TH1F* etaallBkg = (TH1F*)etaspectrumBkg->Clone();
918 
919  TH1F* etapidMass = (TH1F*)etaspectrumMass->Clone();
920  TH1F* etapidSgn = (TH1F*)etaspectrumSgn->Clone();
921  TH1F* etapidBkg = (TH1F*)etaspectrumBkg->Clone();
922 
923  fOutputAll->Add(etaallMass);
924  fOutputAll->Add(etaallSgn);
925  fOutputAll->Add(etaallBkg);
926  fAllhist[(fNPtBins+2)*ketaMass]=etaallMass;
927  fAllhist[(fNPtBins+2)*ketaSgn]=etaallSgn;
928  fAllhist[(fNPtBins+2)*ketaBkg]=etaallBkg;
929 
930  fOutputPID->Add(etapidMass);
931  fOutputPID->Add(etapidSgn);
932  fOutputPID->Add(etapidBkg);
933  fPIDhist[(fNPtBins+2)*ketaMass]=etapidMass;
934  fPIDhist[(fNPtBins+2)*ketaSgn]=etapidSgn;
935  fPIDhist[(fNPtBins+2)*ketaBkg]=etapidBkg;
936 
937  if (fDoDStarVsY){
938  TH3F* deltamassVsyVsPtPID = new TH3F("deltamassVsyVsPt", "delta mass Vs y Vs pT; #DeltaM [GeV/c^{2}]; y; p_{T} [GeV/c]", 700,0.13,0.2, 40, -1, 1, 36, 0., 36.);
939  fOutputPID->Add(deltamassVsyVsPtPID);
940  }
941  return;
942 }
943 //________________________________________________________________________
945  //
947  //
948 
949  if(!isSel) return;
950 
951  // D0 window
952  Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
953  Double_t invmassD0 = part->InvMassD0();
954 
955 
956  Int_t ptbin=cuts->PtBin(part->Pt());
957  Double_t pt = part->Pt();
958  Double_t eta = part->Eta();
959 
960  Double_t invmassDelta = part->DeltaInvMass();
961  Double_t invmassDstar = part->InvMassDstarKpipi();
962 
963  TString fillthis="";
964  Bool_t massInRange=kFALSE;
965 
966  Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
967 
968  // delta M(Kpipi)-M(Kpi)
969  if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<fPeakWindow) massInRange=kTRUE;
970 
971  if(fUseMCInfo) {
972  if(isDStar==1) {
973  histlist[ptbin+1+((fNPtBins+2)*kDzSgn)]->Fill(invmassD0);
974  histlist[(fNPtBins+2)*kDzSgn]->Fill(invmassD0);
975  histlist[ptbin+1+((fNPtBins+2)*kDstarSgn)]->Fill(invmassDstar);
976  histlist[(fNPtBins+2)*kDstarSgn]->Fill(invmassDstar);
977  histlist[ptbin+1+((fNPtBins+2)*kDeltaSgn)]->Fill(invmassDelta);
978  histlist[(fNPtBins+2)*kDeltaSgn]->Fill(invmassDelta);
979  if (massInRange) {
980  histlist[(fNPtBins+2)*kptSgn]->Fill(pt);
981  histlist[(fNPtBins+2)*ketaSgn]->Fill(eta);
982  }
983  }
984  else {//background
985  histlist[ptbin+1+((fNPtBins+2)*kDzBkg)]->Fill(invmassD0);
986  histlist[(fNPtBins+2)*kDzBkg]->Fill(invmassD0);
987  histlist[ptbin+1+((fNPtBins+2)*kDstarBkg)]->Fill(invmassDstar);
988  histlist[(fNPtBins+2)*kDstarBkg]->Fill(invmassDstar);
989  histlist[ptbin+1+((fNPtBins+2)*kDeltaBkg)]->Fill(invmassDelta);
990  histlist[(fNPtBins+2)*kDeltaBkg]->Fill(invmassDelta);
991  if (massInRange) {
992  histlist[(fNPtBins+2)*kptBkg]->Fill(pt);
993  histlist[(fNPtBins+2)*ketaBkg]->Fill(eta);
994  }
995  }
996  }
997  //no MC info, just cut selection
998  histlist[ptbin+1+((fNPtBins+2)*kDzMass)]->Fill(invmassD0);
999  histlist[(fNPtBins+2)*kDzMass]->Fill(invmassD0);
1000  histlist[ptbin+1+((fNPtBins+2)*kDstarMass)]->Fill(invmassDstar);
1001  histlist[(fNPtBins+2)*kDstarMass]->Fill(invmassDstar);
1002  histlist[ptbin+1+((fNPtBins+2)*kDeltaMass)]->Fill(invmassDelta);
1003  histlist[(fNPtBins+2)*kDeltaMass]->Fill(invmassDelta);
1004 
1005  if (massInRange) {
1006  histlist[(fNPtBins+2)*kptMass]->Fill(pt);
1007  histlist[(fNPtBins+2)*ketaMass]->Fill(eta);
1008  }
1009 
1010  return;
1011 }
1012 //______________________________ side band background for D*___________________________________
1014 
1018 
1019  if(!isSel) return;
1020 
1021  Int_t ptbin=cuts->PtBin(part->Pt());
1022 
1023  // select the side bands intervall
1024  Double_t invmassD0 = part->InvMassD0();
1025  if(TMath::Abs(invmassD0-1.865)>4*fD0Window && TMath::Abs(invmassD0-1.865)<8*fD0Window){
1026 
1027  // for pt and eta
1028  Double_t invmassDelta = part->DeltaInvMass();
1029 
1030  histlist[ptbin+1+((fNPtBins+2)*kSideBandMass)]->Fill(invmassDelta);
1031  histlist[(fNPtBins+2)*kSideBandMass]->Fill(invmassDelta);
1032 
1033 
1034  }
1035 }
1036 //________________________________________________________________________________________________________________
1038  //
1040  //
1041  Int_t ptbin=cuts->PtBin(part->Pt());
1042 
1043  Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1044  Double_t invmassD0 = part->InvMassD0();
1045  if (TMath::Abs(invmassD0-mPDGD0)>fD0Window) return;
1046 
1047  AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)part->Get2Prong();
1048 
1049  Int_t okDzWrongSign;
1050  Double_t wrongMassD0=0.;
1051 
1052  Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kCandidate); //selected
1053  if (!isSelected){
1054  return;
1055  }
1056 
1057  okDzWrongSign = 1;
1058 
1059  //if is D*+ than assume D0bar
1060  if(part->Charge()>0 && (isSelected ==1)) {
1061  okDzWrongSign = 0;
1062  }
1063 
1064  // assign the wrong mass in case the cuts return both D0 and D0bar
1065  if(part->Charge()>0 && (isSelected ==3)) {
1066  okDzWrongSign = 0;
1067  }
1068 
1069  //wrong D0 inv mass
1070  if(okDzWrongSign!=0){
1071  wrongMassD0 = theD0particle->InvMassD0();
1072  }else if(okDzWrongSign==0){
1073  wrongMassD0 = theD0particle->InvMassD0bar();
1074  }
1075 
1076  if(TMath::Abs(wrongMassD0-1.865)<fD0Window){
1077 
1078  // wrong D* inv mass
1079  Double_t e[3];
1080  if (part->Charge()>0){
1081  e[0]=theD0particle->EProng(0,321);
1082  e[1]=theD0particle->EProng(1,211);
1083  }else{
1084  e[0]=theD0particle->EProng(0,211);
1085  e[1]=theD0particle->EProng(1,321);
1086  }
1087  e[2]=part->EProng(0,211);
1088 
1089  Double_t esum = e[0]+e[1]+e[2];
1090  Double_t pds = part->P();
1091 
1092  Double_t wrongMassDstar = TMath::Sqrt(esum*esum-pds*pds);
1093 
1094  TString fillthis="";
1095  fillthis="histWrongSignMass_";
1096  fillthis+=ptbin;
1097  ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);
1098  fillthis="histWrongSignMass";
1099  ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);
1100 
1101  }
1102 }
1103 
1104 //-------------------------------------------------------------------------------
1105 Int_t AliAnalysisTaskSEDStarSpectra::CheckOrigin(TClonesArray* arrayMC, const AliAODMCParticle *mcPartCandidate) const {
1106  //
1107  // checking whether the mother of the particles come from a charm or a bottom quark
1108  //
1109 
1110  Int_t pdgGranma = 0;
1111  Int_t mother = 0;
1112  mother = mcPartCandidate->GetMother();
1113  Int_t istep = 0;
1114  Int_t abspdgGranma =0;
1115  Bool_t isFromB=kFALSE;
1116  while (mother >0 ){
1117  istep++;
1118  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1119  if (mcGranma){
1120  pdgGranma = mcGranma->GetPdgCode();
1121  abspdgGranma = TMath::Abs(pdgGranma);
1122  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1123  isFromB=kTRUE;
1124  }
1125  mother = mcGranma->GetMother();
1126  }else{
1127  AliError("Failed casting the mother particle!");
1128  break;
1129  }
1130  }
1131 
1132  if(isFromB) return 5;
1133  else return 4;
1134 }
1135 //-------------------------------------------------------------------------------------
1136 Float_t AliAnalysisTaskSEDStarSpectra::GetTrueImpactParameterD0(const AliAODMCHeader *mcHeader, TClonesArray* arrayMC, const AliAODMCParticle *partDp) const {
1138 
1139  Double_t vtxTrue[3];
1140  mcHeader->GetVertex(vtxTrue);
1141  Double_t origD[3];
1142  partDp->XvYvZv(origD);
1143  Short_t charge=partDp->Charge();
1144  Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
1145  Int_t labelFirstDau = partDp->GetDaughter(0);
1146 
1147  Int_t nDau=partDp->GetNDaughters();
1148 
1149  Int_t theDau=0;
1150  if(nDau==2){
1151  for(Int_t iDau=0; iDau<2; iDau++){
1152  Int_t ind = labelFirstDau+iDau;
1153  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1154  if(!part){
1155  AliError("Daughter particle not found in MC array");
1156  return 99999.;
1157  }
1158  Int_t pdgCode=TMath::Abs(part->GetPdgCode());
1159  if(pdgCode==211 || pdgCode==321){
1160  pXdauTrue[theDau]=part->Px();
1161  pYdauTrue[theDau]=part->Py();
1162  pZdauTrue[theDau]=part->Pz();
1163  ++theDau;
1164  }
1165  }
1166  }
1167  if(theDau!=2){
1168  AliError("Wrong number of decay prongs");
1169  return 99999.;
1170  }
1171 
1172  Double_t d0dummy[3]={0.,0.,0.};
1173  AliAODRecoDecayHF aodD0MC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
1174  return aodD0MC.ImpParXY();
1175 
1176 }
1177 //______________________________________________________-
1180 
1181  Int_t nbins[3]={400,200,fNImpParBins};
1182  Double_t xmin[3]={1.75,0.,fLowerImpPar};
1183  Double_t xmax[3]={1.98,20.,fHigherImpPar};
1184 
1185  fHistMassPtImpParTCDs[0]=new THnSparseF("hMassPtImpParAll",
1186  "Mass vs. pt vs.imppar - All",
1187  3,nbins,xmin,xmax);
1188  fHistMassPtImpParTCDs[1]=new THnSparseF("hMassPtImpParPrompt",
1189  "Mass vs. pt vs.imppar - promptD",
1190  3,nbins,xmin,xmax);
1191  fHistMassPtImpParTCDs[2]=new THnSparseF("hMassPtImpParBfeed",
1192  "Mass vs. pt vs.imppar - DfromB",
1193  3,nbins,xmin,xmax);
1194  fHistMassPtImpParTCDs[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1195  "Mass vs. pt vs.true imppar -DfromB",
1196  3,nbins,xmin,xmax);
1197  fHistMassPtImpParTCDs[4]=new THnSparseF("hMassPtImpParBkg",
1198  "Mass vs. pt vs.imppar - backgr.",
1199  3,nbins,xmin,xmax);
1200 
1201  for(Int_t i=0; i<5;i++){
1202  fOutput->Add(fHistMassPtImpParTCDs[i]);
1203  }
1204 }
1205 
Int_t charge
TH1F * fCEvents
Cuts - sent to output slot 3.
double Double_t
Definition: External.C:58
Definition: External.C:260
void StoreCandidates(AliVEvent *, Int_t nCand=0, Bool_t flagFilter=kTRUE)
Definition: External.C:236
Double_t DeltaInvMass() const
Float_t fHigherImpPar
lower limit in impact parameter (um)
Double_t YDstar() const
Int_t CheckOrigin(TClonesArray *arrayMC, const AliAODMCParticle *mcPartCandidate) const
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
Float_t fLowerImpPar
nunber of bins in impact parameter histos
AliRDHFCutsDStartoKpipi * fCuts
n sigma for kaon PID
Double_t mass
Double_t ImpParXY() const
char Char_t
Definition: External.C:18
static Int_t CheckMatchingAODdeltaAODevents()
Int_t GetWhyRejection() const
Definition: AliRDHFCuts.h:315
Bool_t fDoDStarVsY
higher limit in impact parameter (um)
Float_t GetTrueImpactParameterD0(const AliAODMCHeader *mcHeader, TClonesArray *arrayMC, const AliAODMCParticle *partDp) const
AliNormalizationCounter * fCounter
!Counter for normalization slot 4
Bool_t FillRecoCasc(AliVEvent *event, AliAODRecoCascadeHF *rc, Bool_t isDStar, Bool_t recoSecVtx=kFALSE)
virtual void UserCreateOutputObjects()
Implementation of interface methods.
int Int_t
Definition: External.C:63
void FillSpectrum(AliAODRecoCascadeHF *part, Int_t isDStar, AliRDHFCutsDStartoKpipi *cuts, Int_t isSel, TList *listout, TH1F **histlist)
histos
float Float_t
Definition: External.C:68
virtual void UserExec(Option_t *option)
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel, AliAODEvent *aod)
virtual void Terminate(Option_t *option)
void SideBandBackground(AliAODRecoCascadeHF *part, AliRDHFCutsDStartoKpipi *cuts, Int_t isSel, TList *listout, TH1F **histlist)
Background simulation.
Bool_t fDoImpParDstar
-1: no protection, 0: check AOD/dAOD nEvents only, 1: check AOD/dAOD nEvents + TProcessID ...
Double_t fPeakWindow
select width on D0Mass
short Short_t
Definition: External.C:23
void WrongSignForDStar(AliAODRecoCascadeHF *part, AliRDHFCutsDStartoKpipi *cuts, TList *listout)
Bool_t IsEventSelected(AliVEvent *event)
void StoreEvent(AliVEvent *, AliRDHFCuts *, Bool_t mc=kFALSE, Int_t multiplicity=-9999, Double_t spherocity=-99.)
void SetUsePID(Bool_t flag=kTRUE)
Definition: AliRDHFCuts.h:209
virtual Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y) const
THnSparseF * fHistMassPtImpParTCDs[5]
flag to enable D* vs y
Double_t InvMassD0() const
const char Option_t
Definition: External.C:48
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:249
Double_t InvMassDstarKpipi() const
const Int_t nbins
bool Bool_t
Definition: External.C:53
AliAODRecoDecayHF2Prong * Get2Prong() const
Int_t PtBin(Double_t pt) const
Bool_t fUseMCInfo
select width on DstarMass