AliPhysics  5a28df1 (5a28df1)
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 search window bin by bin - useful to calculate side band bkg
358  if (0<=Dstarpt && Dstarpt<0.5){
359  if(fAnalysis==1){
360  fD0Window=0.035;
361  fPeakWindow=0.03;
362  }else{
363  fD0Window=0.020;
364  fPeakWindow=0.0018;
365  }
366  }
367  if (0.5<=Dstarpt && Dstarpt<1.0){
368  if(fAnalysis==1){
369  fD0Window=0.035;
370  fPeakWindow=0.03;
371  }else{
372  fD0Window=0.020;
373  fPeakWindow=0.0018;
374  }
375  }
376  if (1.0<=Dstarpt && Dstarpt<2.0){
377  if(fAnalysis==1){
378  fD0Window=0.035;
379  fPeakWindow=0.03;
380  }else{
381  fD0Window=0.020;
382  fPeakWindow=0.0018;
383  }
384  }
385  if (2.0<=Dstarpt && Dstarpt<3.0){
386  if(fAnalysis==1){
387  fD0Window=0.035;
388  fPeakWindow=0.03;
389  }else{
390  fD0Window=0.022;
391  fPeakWindow=0.0016;
392  }
393  }
394  if (3.0<=Dstarpt && Dstarpt<4.0){
395  if(fAnalysis==1){
396  fD0Window=0.035;
397  fPeakWindow=0.03;
398  }else{
399  fD0Window=0.026;
400  fPeakWindow=0.0014;
401  }
402  }
403  if (4.0<=Dstarpt && Dstarpt<5.0){
404  if(fAnalysis==1){
405  fD0Window=0.045;
406  fPeakWindow=0.03;
407  }else{
408  fD0Window=0.026;
409  fPeakWindow=0.0014;
410  }
411  }
412  if (5.0<=Dstarpt && Dstarpt<6.0){
413  if(fAnalysis==1){
414  fD0Window=0.045;
415  fPeakWindow=0.03;
416  }else{
417  fD0Window=0.026;
418  fPeakWindow=0.006;
419  }
420  }
421  if (6.0<=Dstarpt && Dstarpt<7.0){
422  if(fAnalysis==1){
423  fD0Window=0.055;
424  fPeakWindow=0.03;
425  }else{
426  fD0Window=0.026;
427  fPeakWindow=0.006;
428  }
429  }
430  if (Dstarpt>=7.0){
431  if(fAnalysis==1){
432  fD0Window=0.074;
433  fPeakWindow=0.03;
434  }else{
435  fD0Window=0.026;
436  fPeakWindow=0.006;
437  }
438  }
439 
440  nSelectedProd++;
441  nSelectedAna++;
442 
443  // check that we are close to signal in the DeltaM - here to save time for PbPb
444  Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
445  Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
446  Double_t invmassDelta = dstarD0pi->DeltaInvMass();
447 
448  if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))>fPeakWindow) continue;
449  Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate,aodEvent); //selected
450 
451  // after cuts
452  if(fDoImpParDstar && isSelected){
453  fHistMassPtImpParTCDs[0]->Fill(arrayForSparse);
454  if(isPrimary) fHistMassPtImpParTCDs[1]->Fill(arrayForSparse);
455  else{
456  fHistMassPtImpParTCDs[2]->Fill(arrayForSparse);
457  fHistMassPtImpParTCDs[3]->Fill(arrayForSparseTrue);
458  }
459  }
460 
461  if (fDoDStarVsY && isSelected){
462  ((TH3F*) (fOutputPID->FindObject("deltamassVsyVsPt")))->Fill(dstarD0pi->DeltaInvMass(),dstarD0pi->YDstar(),dstarD0pi->Pt() );
463  }
464 
465 
466  // fill PID
467  FillSpectrum(dstarD0pi,isDStar,fCuts,isSelected,fOutputPID, fPIDhist);
468  SideBandBackground(dstarD0pi,fCuts,isSelected, fOutputPID, fPIDhist);
469  //WrongSignForDStar(dstarD0pi,fCuts,fOutputPID);
470 
471  //swich off the PID selection
472  fCuts->SetUsePID(kFALSE);
473  Int_t isSelectedNoPID=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate, aodEvent); //selected
474  fCuts->SetUsePID(kTRUE);
475 
476  FillSpectrum(dstarD0pi,isDStar,fCuts,isSelectedNoPID,fOutputAll, fAllhist);
477  // SideBandBackground(dstarD0pi,fCuts,isSelectedNoPID, fOutputAll);
478 
479  // rare D search ------
480  if(fDoSearch){
481  TLorentzVector lorentzTrack1(0,0,0,0); // lorentz 4 vector
482  TLorentzVector lorentzTrack2(0,0,0,0); // lorentz 4 vector
483 
484  for (Int_t i=0; i<aodEvent->GetNumberOfTracks(); i++){
485 
486  AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(i));
487  if(!aodTrack) AliFatal("Not a standard AOD");
488 
489  if(dstarD0pi->Charge() == aodTrack->Charge()) continue;
490  if((!(aodTrack->GetStatus()&AliESDtrack::kITSrefit)|| (!(aodTrack->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
491  if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))>0.02) continue;
492 
493  //build the D1 mass
494  Double_t mass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
495 
496  lorentzTrack1.SetPxPyPzE( dstarD0pi->Px(),dstarD0pi->Py(), dstarD0pi->Pz(), dstarD0pi->E(413) );
497  lorentzTrack2.SetPxPyPzE( aodTrack->Px(),aodTrack->Py(), aodTrack->Pz(),aodTrack->E(mass) );
498 
499  //D1 mass
500  Double_t d1mass = ((lorentzTrack1+lorentzTrack2).M());
501  //mass difference - at 0.4117 and 0.4566
502  fDeltaMassD1->Fill(d1mass-dstarD0pi->InvMassDstarKpipi());
503  }
504  }
505 
506  if(isDStar == 1) {
507  fTrueDiff2->Fill(dstarD0pi->Pt(),dstarD0pi->DeltaInvMass());
508  }
509 
510  }
511 
512  fCounter->StoreCandidates(aodEvent,nSelectedProd,kTRUE);
513  fCounter->StoreCandidates(aodEvent,nSelectedAna,kFALSE);
514 
515  delete vHF;
516 
517  AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
518 
519  PostData(1,fOutput);
520  PostData(2,fOutputAll);
521  PostData(3,fOutputPID);
522  PostData(5,fCounter);
523 
524 }
525 //________________________________________ terminate ___________________________
527 {
531 
532  //Info("Terminate","");
533  AliAnalysisTaskSE::Terminate();
534 
535  fOutput = dynamic_cast<TList*> (GetOutputData(1));
536  if (!fOutput) {
537  printf("ERROR: fOutput not available\n");
538  return;
539  }
540 
541  fCEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fCEvents"));
542  fDeltaMassD1 = dynamic_cast<TH1F*>(fOutput->FindObject("fDeltaMassD1"));
543  fTrueDiff2 = dynamic_cast<TH2F*>(fOutput->FindObject("fTrueDiff2"));
544 
545  fOutputAll = dynamic_cast<TList*> (GetOutputData(1));
546  if (!fOutputAll) {
547  printf("ERROR: fOutputAll not available\n");
548  return;
549  }
550  fOutputPID = dynamic_cast<TList*> (GetOutputData(2));
551  if (!fOutputPID) {
552  printf("ERROR: fOutputPID not available\n");
553  return;
554  }
555  return;
556 }
557 //___________________________________________________________________________
560  Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
561 
562  //slot #1
563  //OpenFile(1);
564  fOutput = new TList();
565  fOutput->SetOwner();
566  fOutput->SetName("chist0");
567 
568  fOutputAll = new TList();
569  fOutputAll->SetOwner();
570  fOutputAll->SetName("listAll");
571 
572  fOutputPID = new TList();
573  fOutputPID->SetOwner();
574  fOutputPID->SetName("listPID");
575 
576  // define histograms
578 
579  //Counter for Normalization
580  fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));
581  fCounter->Init();
582 
584 
585  PostData(1,fOutput);
586  PostData(2,fOutputAll);
587  PostData(3,fOutputPID);
588 
589  return;
590 }
591 //___________________________________ hiostograms _______________________________________
594 
595  fCEvents = new TH1F("fCEvents","counter",13,0,13);
596  fCEvents->SetStats(kTRUE);
597  fCEvents->GetXaxis()->SetTitle("1");
598  fCEvents->GetYaxis()->SetTitle("counts");
599  fCEvents->GetXaxis()->SetBinLabel(2,"no. of events");
600  fCEvents->GetXaxis()->SetBinLabel(3,"good prim vtx and B field");
601  fCEvents->GetXaxis()->SetBinLabel(4,"no event selected");
602  fCEvents->GetXaxis()->SetBinLabel(5,"no vtx contributors");
603  fCEvents->GetXaxis()->SetBinLabel(6,"trigger for PbPb");
604  fCEvents->GetXaxis()->SetBinLabel(7,"no z vtx");
605  fCEvents->GetXaxis()->SetBinLabel(12,"no. of D0 fail to be rec");
606  fOutput->Add(fCEvents);
607 
608  fTrueDiff2 = new TH2F("DiffDstar_pt","True Reco diff vs pt",200,0,15,900,0,0.3);
609  fOutput->Add(fTrueDiff2);
610 
611  fDeltaMassD1 = new TH1F("DeltaMassD1","delta mass d1",600,0,0.8);
612  fOutput->Add(fDeltaMassD1);
613  //temp a
614 
615  fAllhist = new TH1F*[(fNPtBins+2)*18];
616  fPIDhist = new TH1F*[(fNPtBins+2)*18];
617 
618  TString nameMass=" ", nameSgn=" ", nameBkg=" ";
619 
620  for(Int_t i=-2;i<fNPtBins;i++){
621  nameMass="histDeltaMass_";
622  nameMass+=i+1;
623  nameSgn="histDeltaSgn_";
624  nameSgn+=i+1;
625  nameBkg="histDeltaBkg_";
626  nameBkg+=i+1;
627 
628  if (i==-2) {
629  nameMass="histDeltaMass";
630  nameSgn="histDeltaSgn";
631  nameBkg="histDeltaBkg";
632  }
633 
634  TH1F* spectrumMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} invariant mass; #DeltaM [GeV/c^{2}]; Entries",700,0.13,0.2);
635  TH1F* spectrumSgn = new TH1F(nameSgn.Data(), "D^{*}-D^{0} Signal invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",700,0.13,0.2);
636  TH1F* spectrumBkg = new TH1F(nameBkg.Data(), "D^{*}-D^{0} Background invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",700,0.13,0.2);
637 
638  nameMass="histD0Mass_";
639  nameMass+=i+1;
640  nameSgn="histD0Sgn_";
641  nameSgn+=i+1;
642  nameBkg="histD0Bkg_";
643  nameBkg+=i+1;
644 
645  if (i==-2) {
646  nameMass="histD0Mass";
647  nameSgn="histD0Sgn";
648  nameBkg="histD0Bkg";
649  }
650 
651  TH1F* spectrumD0Mass = new TH1F(nameMass.Data(),"D^{0} invariant mass; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
652  TH1F* spectrumD0Sgn = new TH1F(nameSgn.Data(), "D^{0} Signal invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
653  TH1F* spectrumD0Bkg = new TH1F(nameBkg.Data(), "D^{0} Background invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
654 
655  nameMass="histDstarMass_";
656  nameMass+=i+1;
657  nameSgn="histDstarSgn_";
658  nameSgn+=i+1;
659  nameBkg="histDstarBkg_";
660  nameBkg+=i+1;
661 
662  if (i==-2) {
663  nameMass="histDstarMass";
664  nameSgn="histDstarSgn";
665  nameBkg="histDstarBkg";
666  }
667 
668  TH1F* spectrumDstarMass = new TH1F(nameMass.Data(),"D^{*} invariant mass; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
669  TH1F* spectrumDstarSgn = new TH1F(nameSgn.Data(), "D^{*} Signal invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
670  TH1F* spectrumDstarBkg = new TH1F(nameBkg.Data(), "D^{*} Background invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
671 
672  nameMass="histSideBandMass_";
673  nameMass+=i+1;
674  if (i==-2) {
675  nameMass="histSideBandMass";
676  }
677 
678  TH1F* spectrumSideBandMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} sideband mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
679 
680  nameMass="histWrongSignMass_";
681  nameMass+=i+1;
682  if (i==-2) {
683  nameMass="histWrongSignMass";
684  }
685 
686  TH1F* spectrumWrongSignMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} wrongsign mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
687 
688 
689  spectrumMass->Sumw2();
690  spectrumSgn->Sumw2();
691  spectrumBkg->Sumw2();
692 
693  spectrumMass->SetLineColor(6);
694  spectrumSgn->SetLineColor(2);
695  spectrumBkg->SetLineColor(4);
696 
697  spectrumMass->SetMarkerStyle(20);
698  spectrumSgn->SetMarkerStyle(20);
699  spectrumBkg->SetMarkerStyle(20);
700  spectrumMass->SetMarkerSize(0.6);
701  spectrumSgn->SetMarkerSize(0.6);
702  spectrumBkg->SetMarkerSize(0.6);
703  spectrumMass->SetMarkerColor(6);
704  spectrumSgn->SetMarkerColor(2);
705  spectrumBkg->SetMarkerColor(4);
706 
707  spectrumD0Mass->Sumw2();
708  spectrumD0Sgn->Sumw2();
709  spectrumD0Bkg->Sumw2();
710 
711  spectrumD0Mass->SetLineColor(6);
712  spectrumD0Sgn->SetLineColor(2);
713  spectrumD0Bkg->SetLineColor(4);
714 
715  spectrumD0Mass->SetMarkerStyle(20);
716  spectrumD0Sgn->SetMarkerStyle(20);
717  spectrumD0Bkg->SetMarkerStyle(20);
718  spectrumD0Mass->SetMarkerSize(0.6);
719  spectrumD0Sgn->SetMarkerSize(0.6);
720  spectrumD0Bkg->SetMarkerSize(0.6);
721  spectrumD0Mass->SetMarkerColor(6);
722  spectrumD0Sgn->SetMarkerColor(2);
723  spectrumD0Bkg->SetMarkerColor(4);
724 
725  spectrumDstarMass->Sumw2();
726  spectrumDstarSgn->Sumw2();
727  spectrumDstarBkg->Sumw2();
728 
729  spectrumDstarMass->SetLineColor(6);
730  spectrumDstarSgn->SetLineColor(2);
731  spectrumDstarBkg->SetLineColor(4);
732 
733  spectrumDstarMass->SetMarkerStyle(20);
734  spectrumDstarSgn->SetMarkerStyle(20);
735  spectrumDstarBkg->SetMarkerStyle(20);
736  spectrumDstarMass->SetMarkerSize(0.6);
737  spectrumDstarSgn->SetMarkerSize(0.6);
738  spectrumDstarBkg->SetMarkerSize(0.6);
739  spectrumDstarMass->SetMarkerColor(6);
740  spectrumDstarSgn->SetMarkerColor(2);
741  spectrumDstarBkg->SetMarkerColor(4);
742 
743  spectrumSideBandMass->Sumw2();
744  spectrumSideBandMass->SetLineColor(4);
745  spectrumSideBandMass->SetMarkerStyle(20);
746  spectrumSideBandMass->SetMarkerSize(0.6);
747  spectrumSideBandMass->SetMarkerColor(4);
748 
749  spectrumWrongSignMass->Sumw2();
750  spectrumWrongSignMass->SetLineColor(4);
751  spectrumWrongSignMass->SetMarkerStyle(20);
752  spectrumWrongSignMass->SetMarkerSize(0.6);
753  spectrumWrongSignMass->SetMarkerColor(4);
754 
755  TH1F* allMass = (TH1F*)spectrumMass->Clone();
756  TH1F* allSgn = (TH1F*)spectrumSgn->Clone();
757  TH1F* allBkg = (TH1F*)spectrumBkg->Clone();
758 
759  TH1F* pidMass = (TH1F*)spectrumMass->Clone();
760  TH1F* pidSgn = (TH1F*)spectrumSgn->Clone();
761  TH1F* pidBkg = (TH1F*)spectrumBkg->Clone();
762 
763  fOutputAll->Add(allMass);
764  fOutputAll->Add(allSgn);
765  fOutputAll->Add(allBkg);
766  fAllhist[i+2+((fNPtBins+2)*kDeltaMass)]=allMass;
767  fAllhist[i+2+((fNPtBins+2)*kDeltaSgn)]=allSgn;
768  fAllhist[i+2+((fNPtBins+2)*kDeltaBkg)]=allBkg;
769 
770  fOutputPID->Add(pidMass);
771  fOutputPID->Add(pidSgn);
772  fOutputPID->Add(pidBkg);
773  fPIDhist[i+2+((fNPtBins+2)*kDeltaMass)]=pidMass;
774  fPIDhist[i+2+((fNPtBins+2)*kDeltaSgn)]=pidSgn;
775  fPIDhist[i+2+((fNPtBins+2)*kDeltaBkg)]=pidBkg;
776 
777  TH1F* allD0Mass = (TH1F*)spectrumD0Mass->Clone();
778  TH1F* allD0Sgn = (TH1F*)spectrumD0Sgn->Clone();
779  TH1F* allD0Bkg = (TH1F*)spectrumD0Bkg->Clone();
780 
781  TH1F* pidD0Mass = (TH1F*)spectrumD0Mass->Clone();
782  TH1F* pidD0Sgn = (TH1F*)spectrumD0Sgn->Clone();
783  TH1F* pidD0Bkg = (TH1F*)spectrumD0Bkg->Clone();
784 
785  fOutputAll->Add(allD0Mass);
786  fOutputAll->Add(allD0Sgn);
787  fOutputAll->Add(allD0Bkg);
788  fAllhist[i+2+((fNPtBins+2)*kDzMass)]=allD0Mass;
789  fAllhist[i+2+((fNPtBins+2)*kDzSgn)]=allD0Sgn;
790  fAllhist[i+2+((fNPtBins+2)*kDzBkg)]=allD0Bkg;
791 
792  fOutputPID->Add(pidD0Mass);
793  fOutputPID->Add(pidD0Sgn);
794  fOutputPID->Add(pidD0Bkg);
795  fPIDhist[i+2+((fNPtBins+2)*kDzMass)]=pidD0Mass;
796  fPIDhist[i+2+((fNPtBins+2)*kDzSgn)]=pidD0Sgn;
797  fPIDhist[i+2+((fNPtBins+2)*kDzBkg)]=pidD0Bkg;
798 
799  TH1F* allDstarMass = (TH1F*)spectrumDstarMass->Clone();
800  TH1F* allDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
801  TH1F* allDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
802 
803  TH1F* pidDstarMass = (TH1F*)spectrumDstarMass->Clone();
804  TH1F* pidDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
805  TH1F* pidDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
806 
807  fOutputAll->Add(allDstarMass);
808  fOutputAll->Add(allDstarSgn);
809  fOutputAll->Add(allDstarBkg);
810  fAllhist[i+2+((fNPtBins+2)*kDstarMass)]=allDstarMass;
811  fAllhist[i+2+((fNPtBins+2)*kDstarSgn)]=allDstarSgn;
812  fAllhist[i+2+((fNPtBins+2)*kDstarBkg)]=allDstarBkg;
813 
814  fOutputPID->Add(pidDstarMass);
815  fOutputPID->Add(pidDstarSgn);
816  fOutputPID->Add(pidDstarBkg);
817  fPIDhist[i+2+((fNPtBins+2)*kDstarMass)]=pidDstarMass;
818  fPIDhist[i+2+((fNPtBins+2)*kDstarSgn)]=pidDstarSgn;
819  fPIDhist[i+2+((fNPtBins+2)*kDstarBkg)]=pidDstarBkg;
820 
821  TH1F* allSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
822  TH1F* pidSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
823 
824  fOutputAll->Add(allSideBandMass);
825  fOutputPID->Add(pidSideBandMass);
826  fAllhist[i+2+((fNPtBins+2)*kSideBandMass)]=allSideBandMass;
827  fPIDhist[i+2+((fNPtBins+2)*kSideBandMass)]=pidSideBandMass;
828 
829  TH1F* allWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
830  TH1F* pidWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
831 
832  fOutputAll->Add(allWrongSignMass);
833  fOutputPID->Add(pidWrongSignMass);
834  fAllhist[i+2+((fNPtBins+2)*kWrongSignMass)]=allWrongSignMass;
835  fPIDhist[i+2+((fNPtBins+2)*kWrongSignMass)]=pidWrongSignMass;
836 
837  }
838 
839  // pt spectra
840  nameMass="ptMass";
841  nameSgn="ptSgn";
842  nameBkg="ptBkg";
843 
844  TH1F* ptspectrumMass = new TH1F(nameMass.Data(),"D^{*} p_{T}; p_{T} [GeV]; Entries",400,0,50);
845  TH1F* ptspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal p_{T} - MC; p_{T} [GeV]; Entries",400,0,50);
846  TH1F* ptspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background p_{T} - MC; p_{T} [GeV]; Entries",400,0,50);
847 
848  ptspectrumMass->Sumw2();
849  ptspectrumSgn->Sumw2();
850  ptspectrumBkg->Sumw2();
851 
852  ptspectrumMass->SetLineColor(6);
853  ptspectrumSgn->SetLineColor(2);
854  ptspectrumBkg->SetLineColor(4);
855 
856  ptspectrumMass->SetMarkerStyle(20);
857  ptspectrumSgn->SetMarkerStyle(20);
858  ptspectrumBkg->SetMarkerStyle(20);
859  ptspectrumMass->SetMarkerSize(0.6);
860  ptspectrumSgn->SetMarkerSize(0.6);
861  ptspectrumBkg->SetMarkerSize(0.6);
862  ptspectrumMass->SetMarkerColor(6);
863  ptspectrumSgn->SetMarkerColor(2);
864  ptspectrumBkg->SetMarkerColor(4);
865 
866  TH1F* ptallMass = (TH1F*)ptspectrumMass->Clone();
867  TH1F* ptallSgn = (TH1F*)ptspectrumSgn->Clone();
868  TH1F* ptallBkg = (TH1F*)ptspectrumBkg->Clone();
869 
870  TH1F* ptpidMass = (TH1F*)ptspectrumMass->Clone();
871  TH1F* ptpidSgn = (TH1F*)ptspectrumSgn->Clone();
872  TH1F* ptpidBkg = (TH1F*)ptspectrumBkg->Clone();
873 
874  fOutputAll->Add(ptallMass);
875  fOutputAll->Add(ptallSgn);
876  fOutputAll->Add(ptallBkg);
877  fAllhist[((fNPtBins+2)*kptMass)]=ptallMass;
878  fAllhist[((fNPtBins+2)*kptSgn)]=ptallSgn;
879  fAllhist[((fNPtBins+2)*kptBkg)]=ptallBkg;
880 
881  fOutputPID->Add(ptpidMass);
882  fOutputPID->Add(ptpidSgn);
883  fOutputPID->Add(ptpidBkg);
884  fPIDhist[(fNPtBins+2)*kptMass]=ptpidMass;
885  fPIDhist[(fNPtBins+2)*kptSgn]=ptpidSgn;
886  fPIDhist[(fNPtBins+2)*kptBkg]=ptpidBkg;
887  // eta spectra
888  nameMass="etaMass";
889  nameSgn="etaSgn";
890  nameBkg="etaBkg";
891 
892  TH1F* etaspectrumMass = new TH1F(nameMass.Data(),"D^{*} #eta; #eta; Entries",200,-1,1);
893  TH1F* etaspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal #eta - MC; #eta; Entries",200,-1,1);
894  TH1F* etaspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background #eta - MC; #eta; Entries",200,-1,1);
895 
896  etaspectrumMass->Sumw2();
897  etaspectrumSgn->Sumw2();
898  etaspectrumBkg->Sumw2();
899 
900  etaspectrumMass->SetLineColor(6);
901  etaspectrumSgn->SetLineColor(2);
902  etaspectrumBkg->SetLineColor(4);
903 
904  etaspectrumMass->SetMarkerStyle(20);
905  etaspectrumSgn->SetMarkerStyle(20);
906  etaspectrumBkg->SetMarkerStyle(20);
907  etaspectrumMass->SetMarkerSize(0.6);
908  etaspectrumSgn->SetMarkerSize(0.6);
909  etaspectrumBkg->SetMarkerSize(0.6);
910  etaspectrumMass->SetMarkerColor(6);
911  etaspectrumSgn->SetMarkerColor(2);
912  etaspectrumBkg->SetMarkerColor(4);
913 
914  TH1F* etaallMass = (TH1F*)etaspectrumMass->Clone();
915  TH1F* etaallSgn = (TH1F*)etaspectrumSgn->Clone();
916  TH1F* etaallBkg = (TH1F*)etaspectrumBkg->Clone();
917 
918  TH1F* etapidMass = (TH1F*)etaspectrumMass->Clone();
919  TH1F* etapidSgn = (TH1F*)etaspectrumSgn->Clone();
920  TH1F* etapidBkg = (TH1F*)etaspectrumBkg->Clone();
921 
922  fOutputAll->Add(etaallMass);
923  fOutputAll->Add(etaallSgn);
924  fOutputAll->Add(etaallBkg);
925  fAllhist[(fNPtBins+2)*ketaMass]=etaallMass;
926  fAllhist[(fNPtBins+2)*ketaSgn]=etaallSgn;
927  fAllhist[(fNPtBins+2)*ketaBkg]=etaallBkg;
928 
929  fOutputPID->Add(etapidMass);
930  fOutputPID->Add(etapidSgn);
931  fOutputPID->Add(etapidBkg);
932  fPIDhist[(fNPtBins+2)*ketaMass]=etapidMass;
933  fPIDhist[(fNPtBins+2)*ketaSgn]=etapidSgn;
934  fPIDhist[(fNPtBins+2)*ketaBkg]=etapidBkg;
935 
936  if (fDoDStarVsY){
937  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.);
938  fOutputPID->Add(deltamassVsyVsPtPID);
939  }
940  return;
941 }
942 //________________________________________________________________________
944  //
946  //
947 
948  if(!isSel) return;
949 
950  // D0 window
951  Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
952  Double_t invmassD0 = part->InvMassD0();
953  if (TMath::Abs(invmassD0-mPDGD0)>fD0Window) return;
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:313
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
TString part
use mixed event to constrain combinatorial background
Definition: InvMassFit.C:52
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:248
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