AliPhysics  8b695ca (8b695ca)
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  fCEvents->Fill(0);//all events
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(8);
204  return;
205  }
206  fCEvents->Fill(1);
207  }
208 
209  fEvents++;
210 
211  AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
212  TClonesArray *arrayDStartoD0pi=0;
213  TClonesArray *arrayD0toKpi=0;
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  arrayD0toKpi=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
228  }
229  } else {
230  arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
231  arrayD0toKpi=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
232  }
233 
234  // fix for temporary bug in ESDfilter
235  // the AODs with null vertex pointer didn't pass the PhysSel
236  if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
237  fCEvents->Fill(2);
238 
239  fCounter->StoreEvent(aodEvent,fCuts,fUseMCInfo);
240 
241  // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD
242  TString trigclass=aodEvent->GetFiredTriggerClasses();
243  if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD")||trigclass.Contains("C0SMH-B-NOPF-ALL")) fCEvents->Fill(5);
244 
245  if(!fCuts->IsEventSelected(aodEvent)) {
246  if(fCuts->GetWhyRejection()==6) // rejected for Z vertex
247  fCEvents->Fill(6);
248  return;
249  }
250 
251  Bool_t isEvSel=fCuts->IsEventSelected(aodEvent);
252  fCEvents->Fill(3);
253  if(!isEvSel) return;
254 
255  // Load the event
256  // AliInfo(Form("Event %d",fEvents));
257  //if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
258 
259  // counters for efficiencies
260  Int_t icountReco = 0;
261 
262  //D* and D0 prongs needed to MatchToMC method
263  Int_t pdgDgDStartoD0pi[2]={421,211};
264  Int_t pdgDgD0toKpi[2]={321,211};
265 
266  // AOD primary vertex
267  AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
268  if(!vtx1) return;
269  if(vtx1->GetNContributors()<1) return;
270  fCEvents->Fill(4);
271 
272  if (!arrayDStartoD0pi || !arrayD0toKpi){
273  AliInfo("Could not find array of HF vertices, skipping the event");
274  return;
275  }else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
276 
277  Int_t nSelectedAna =0;
278  Int_t nSelectedProd =0;
279 
280  // vHF object is needed to call the method that refills the missing info of the candidates
281  // if they have been deleted in dAOD reconstruction phase
282  // in order to reduce the size of the file
284 
285  // loop over the tracks to search for candidates soft pion
286  for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
287 
288  // D* candidates and D0 from D*
289  AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
290  AliAODRecoDecayHF2Prong *trackD0=(AliAODRecoDecayHF2Prong*)arrayD0toKpi->At(dstarD0pi->GetProngID(1));
291 
292  fCEvents->Fill(10);
293  TObjArray arrTracks(3);
294  for(Int_t ipr=0;ipr<3;ipr++){
295  AliAODTrack *tr;
296  if(ipr == 0) tr=vHF->GetProng(aodEvent,dstarD0pi,ipr); //soft pion
297  else tr=vHF->GetProng(aodEvent,trackD0,ipr-1); //D0 daughters
298  arrTracks.AddAt(tr,ipr);
299  }
300  if(!fCuts->PreSelect(arrTracks)){
301  fCEvents->Fill(13);
302  continue;
303  }
304 
305  Bool_t isDStarCand =kTRUE;
306  if(!(vHF->FillRecoCasc(aodEvent,dstarD0pi,isDStarCand))) {//Fill the data members of the candidate only if they are empty.
307  fCEvents->Fill(12); //monitor how often this fails
308  continue;
309  }
310  if(!dstarD0pi->GetSecondaryVtx()) continue;
311  AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
312  if (!theD0particle) continue;
313 
314  Int_t isDStar = 0;
315  TClonesArray *mcArray = 0;
316  AliAODMCHeader *mcHeader=0;
317 
318  Bool_t isPrimary=kTRUE;
319  Float_t pdgCode=-2;
320  Float_t trueImpParXY=0.;
321 
322  // mc analysis
323  if(fUseMCInfo){
324  //MC array need for maching
325  mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
326  if (!mcArray) {
327  AliError("Could not find Monte-Carlo in AOD");
328  return;
329  }
330  // load MC header
331  mcHeader = (AliAODMCHeader*)aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
332  if(!mcHeader) {
333  printf("AliAnalysisTaskSEDplus::UserExec: MC header branch not found!\n");
334  return;
335  }
336  // find associated MC particle for D* ->D0toKpi
337  Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
338  if(mcLabel>=0){
339 
340  AliAODMCParticle *partDSt = (AliAODMCParticle*)mcArray->At(mcLabel);
341  Int_t checkOrigin = CheckOrigin(mcArray,partDSt);
342  if(checkOrigin==5) isPrimary=kFALSE;
343  AliAODMCParticle *dg0 = (AliAODMCParticle*)mcArray->At(partDSt->GetDaughter(0));
344  // AliAODMCParticle *dg01 = (AliAODMCParticle*)mcArray->At(dg0->GetDaughter(0));
345  pdgCode=TMath::Abs(partDSt->GetPdgCode());
346  if(!isPrimary){
347  trueImpParXY=GetTrueImpactParameterD0(mcHeader,mcArray,dg0)*1000.;
348  }
349  isDStar = 1;
350  }else{
351  pdgCode=-1;
352  }
353  }
354 
355  if(pdgCode==-1) AliDebug(2,"No particle assigned! check\n");
356 
357  Double_t Dstarpt = dstarD0pi->Pt();
358 
359  // quality selction on tracks and region of interest
360  Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
361  if(!isTkSelected) continue;
362 
363  if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
364 
365 
366  //histos for impact par studies - D0!!!
367  Double_t ptCand = dstarD0pi->Get2Prong()->Pt();
368  Double_t invMass=dstarD0pi->InvMassD0();
369  Double_t impparXY=dstarD0pi->Get2Prong()->ImpParXY()*10000.;
370 
371  Double_t arrayForSparse[3]={invMass,ptCand,impparXY};
372  Double_t arrayForSparseTrue[3]={invMass,ptCand,trueImpParXY};
373 
374  // 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
375 
376  if (0<=Dstarpt && Dstarpt<0.5){
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 (0.5<=Dstarpt && Dstarpt<1.0){
386  if(fAnalysis==1){
387  fD0Window=0.035;
388  fPeakWindow=0.03;
389  }else{
390  fD0Window=0.020;
391  fPeakWindow=0.0018;
392  }
393  }
394  if (1.0<=Dstarpt && Dstarpt<2.0){
395  if(fAnalysis==1){
396  fD0Window=0.035;
397  fPeakWindow=0.03;
398  }else{
399  fD0Window=0.020;
400  fPeakWindow=0.0018;
401  }
402  }
403  if (2.0<=Dstarpt && Dstarpt<3.0){
404  if(fAnalysis==1){
405  fD0Window=0.035;
406  fPeakWindow=0.03;
407  }else{
408  fD0Window=0.022;
409  fPeakWindow=0.0016;
410  }
411  }
412  if (3.0<=Dstarpt && Dstarpt<4.0){
413  if(fAnalysis==1){
414  fD0Window=0.035;
415  fPeakWindow=0.03;
416  }else{
417  fD0Window=0.026;
418  fPeakWindow=0.0014;
419  }
420  }
421  if (4.0<=Dstarpt && Dstarpt<5.0){
422  if(fAnalysis==1){
423  fD0Window=0.045;
424  fPeakWindow=0.03;
425  }else{
426  fD0Window=0.026;
427  fPeakWindow=0.0014;
428  }
429  }
430  if (5.0<=Dstarpt && Dstarpt<6.0){
431  if(fAnalysis==1){
432  fD0Window=0.045;
433  fPeakWindow=0.03;
434  }else{
435  fD0Window=0.026;
436  fPeakWindow=0.006;
437  }
438  }
439  if (6.0<=Dstarpt && Dstarpt<7.0){
440  if(fAnalysis==1){
441  fD0Window=0.055;
442  fPeakWindow=0.03;
443  }else{
444  fD0Window=0.026;
445  fPeakWindow=0.006;
446  }
447  }
448  if (Dstarpt>=7.0){
449  if(fAnalysis==1){
450  fD0Window=0.074;
451  fPeakWindow=0.03;
452  }else{
453  fD0Window=0.026;
454  fPeakWindow=0.006;
455  }
456  }
457 
458  nSelectedProd++;
459  nSelectedAna++;
460 
461  // check that we are close to signal in the DeltaM - here to save time for PbPb
462  Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
463  Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
464  Double_t invmassDelta = dstarD0pi->DeltaInvMass();
465 
466  if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))>fPeakWindow) continue;
467  Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate,aodEvent); //selected
468  if(isSelected>0) fCEvents->Fill(11);
469 
470  // after cuts
471  if(fDoImpParDstar && isSelected){
472  fHistMassPtImpParTCDs[0]->Fill(arrayForSparse);
473  if(isPrimary) fHistMassPtImpParTCDs[1]->Fill(arrayForSparse);
474  else{
475  fHistMassPtImpParTCDs[2]->Fill(arrayForSparse);
476  fHistMassPtImpParTCDs[3]->Fill(arrayForSparseTrue);
477  }
478  }
479 
480  if (fDoDStarVsY && isSelected){
481  ((TH3F*) (fOutputPID->FindObject("deltamassVsyVsPt")))->Fill(dstarD0pi->DeltaInvMass(),dstarD0pi->YDstar(),dstarD0pi->Pt() );
482  }
483 
484 
485  // fill PID
486  FillSpectrum(dstarD0pi,isDStar,fCuts,isSelected,fOutputPID, fPIDhist);
487  SideBandBackground(dstarD0pi,fCuts,isSelected, fOutputPID, fPIDhist);
488  //WrongSignForDStar(dstarD0pi,fCuts,fOutputPID);
489 
490  //swich off the PID selection
491  fCuts->SetUsePID(kFALSE);
492  Int_t isSelectedNoPID=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate, aodEvent); //selected
493  fCuts->SetUsePID(kTRUE);
494 
495  FillSpectrum(dstarD0pi,isDStar,fCuts,isSelectedNoPID,fOutputAll, fAllhist);
496  // SideBandBackground(dstarD0pi,fCuts,isSelectedNoPID, fOutputAll);
497 
498  // rare D search ------
499  if(fDoSearch){
500  TLorentzVector lorentzTrack1(0,0,0,0); // lorentz 4 vector
501  TLorentzVector lorentzTrack2(0,0,0,0); // lorentz 4 vector
502 
503  for (Int_t i=0; i<aodEvent->GetNumberOfTracks(); i++){
504 
505  AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(i));
506  if(!aodTrack) AliFatal("Not a standard AOD");
507 
508  if(dstarD0pi->Charge() == aodTrack->Charge()) continue;
509  if((!(aodTrack->GetStatus()&AliESDtrack::kITSrefit)|| (!(aodTrack->GetStatus()&AliESDtrack::kTPCrefit)))) continue;
510  if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))>0.02) continue;
511 
512  //build the D1 mass
513  Double_t mass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
514 
515  lorentzTrack1.SetPxPyPzE( dstarD0pi->Px(),dstarD0pi->Py(), dstarD0pi->Pz(), dstarD0pi->E(413) );
516  lorentzTrack2.SetPxPyPzE( aodTrack->Px(),aodTrack->Py(), aodTrack->Pz(),aodTrack->E(mass) );
517 
518  //D1 mass
519  Double_t d1mass = ((lorentzTrack1+lorentzTrack2).M());
520  //mass difference - at 0.4117 and 0.4566
521  fDeltaMassD1->Fill(d1mass-dstarD0pi->InvMassDstarKpipi());
522  }
523  }
524 
525  if(isDStar == 1) {
526  fTrueDiff2->Fill(dstarD0pi->Pt(),dstarD0pi->DeltaInvMass());
527  }
528 
529  }
530 
531  fCounter->StoreCandidates(aodEvent,nSelectedProd,kTRUE);
532  fCounter->StoreCandidates(aodEvent,nSelectedAna,kFALSE);
533 
534  delete vHF;
535 
536  AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
537 
538  PostData(1,fOutput);
539  PostData(2,fOutputAll);
540  PostData(3,fOutputPID);
541  PostData(5,fCounter);
542 
543 }
544 //________________________________________ terminate ___________________________
546 {
550 
551  //Info("Terminate","");
552  AliAnalysisTaskSE::Terminate();
553 
554  fOutput = dynamic_cast<TList*> (GetOutputData(1));
555  if (!fOutput) {
556  printf("ERROR: fOutput not available\n");
557  return;
558  }
559 
560  fCEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fCEvents"));
561  fDeltaMassD1 = dynamic_cast<TH1F*>(fOutput->FindObject("fDeltaMassD1"));
562  fTrueDiff2 = dynamic_cast<TH2F*>(fOutput->FindObject("fTrueDiff2"));
563 
564  fOutputAll = dynamic_cast<TList*> (GetOutputData(1));
565  if (!fOutputAll) {
566  printf("ERROR: fOutputAll not available\n");
567  return;
568  }
569  fOutputPID = dynamic_cast<TList*> (GetOutputData(2));
570  if (!fOutputPID) {
571  printf("ERROR: fOutputPID not available\n");
572  return;
573  }
574  return;
575 }
576 //___________________________________________________________________________
579  Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
580 
581  //slot #1
582  //OpenFile(1);
583  fOutput = new TList();
584  fOutput->SetOwner();
585  fOutput->SetName("chist0");
586 
587  fOutputAll = new TList();
588  fOutputAll->SetOwner();
589  fOutputAll->SetName("listAll");
590 
591  fOutputPID = new TList();
592  fOutputPID->SetOwner();
593  fOutputPID->SetName("listPID");
594 
595  // define histograms
597 
598  //Counter for Normalization
599  fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));
600  fCounter->Init();
601 
603 
604  PostData(1,fOutput);
605  PostData(2,fOutputAll);
606  PostData(3,fOutputPID);
607 
608  return;
609 }
610 //___________________________________ hiostograms _______________________________________
613 
614  fCEvents = new TH1F("fCEvents","counter",14,0,14);
615  fCEvents->SetStats(kTRUE);
616  fCEvents->GetXaxis()->SetTitle("1");
617  fCEvents->GetYaxis()->SetTitle("counts");
618  fCEvents->GetXaxis()->SetBinLabel(1,"nEventsRead");
619  fCEvents->GetXaxis()->SetBinLabel(2,"nEvents Matched dAOD");
620  fCEvents->GetXaxis()->SetBinLabel(3,"good prim vtx and B field");
621  fCEvents->GetXaxis()->SetBinLabel(4,"no event selected");
622  fCEvents->GetXaxis()->SetBinLabel(5,"no vtx contributors");
623  fCEvents->GetXaxis()->SetBinLabel(6,"trigger for PbPb");
624  fCEvents->GetXaxis()->SetBinLabel(7,"no z vtx");
625  fCEvents->GetXaxis()->SetBinLabel(9,"nEvents Mismatched dAOD");
626  fCEvents->GetXaxis()->SetBinLabel(11, "no. of cascade candidates");
627  fCEvents->GetXaxis()->SetBinLabel(12, "no. of Dstar after selection cuts");
628  fCEvents->GetXaxis()->SetBinLabel(13, "no. of not on-the-fly rec Dstar");
629  fCEvents->GetXaxis()->SetBinLabel(14, "no. of Dstar rejected by preselect"); //toadd
630 
631  fOutput->Add(fCEvents);
632 
633  fTrueDiff2 = new TH2F("DiffDstar_pt","True Reco diff vs pt",200,0,15,900,0,0.3);
634  fOutput->Add(fTrueDiff2);
635 
636  fDeltaMassD1 = new TH1F("DeltaMassD1","delta mass d1",600,0,0.8);
637  fOutput->Add(fDeltaMassD1);
638  //temp a
639 
640  fAllhist = new TH1F*[(fNPtBins+2)*18];
641  fPIDhist = new TH1F*[(fNPtBins+2)*18];
642 
643  TString nameMass=" ", nameSgn=" ", nameBkg=" ";
644 
645  for(Int_t i=-2;i<fNPtBins;i++){
646  nameMass="histDeltaMass_";
647  nameMass+=i+1;
648  nameSgn="histDeltaSgn_";
649  nameSgn+=i+1;
650  nameBkg="histDeltaBkg_";
651  nameBkg+=i+1;
652 
653  if (i==-2) {
654  nameMass="histDeltaMass";
655  nameSgn="histDeltaSgn";
656  nameBkg="histDeltaBkg";
657  }
658 
659  TH1F* spectrumMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} invariant mass; #DeltaM [GeV/c^{2}]; Entries",700,0.13,0.2);
660  TH1F* spectrumSgn = new TH1F(nameSgn.Data(), "D^{*}-D^{0} Signal invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",700,0.13,0.2);
661  TH1F* spectrumBkg = new TH1F(nameBkg.Data(), "D^{*}-D^{0} Background invariant mass - MC; #DeltaM [GeV/c^{2}]; Entries",700,0.13,0.2);
662 
663  nameMass="histD0Mass_";
664  nameMass+=i+1;
665  nameSgn="histD0Sgn_";
666  nameSgn+=i+1;
667  nameBkg="histD0Bkg_";
668  nameBkg+=i+1;
669 
670  if (i==-2) {
671  nameMass="histD0Mass";
672  nameSgn="histD0Sgn";
673  nameBkg="histD0Bkg";
674  }
675 
676  TH1F* spectrumD0Mass = new TH1F(nameMass.Data(),"D^{0} invariant mass; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
677  TH1F* spectrumD0Sgn = new TH1F(nameSgn.Data(), "D^{0} Signal invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
678  TH1F* spectrumD0Bkg = new TH1F(nameBkg.Data(), "D^{0} Background invariant mass - MC; M(D^{0}) [GeV/c^{2}]; Entries",200,1.75,1.95);
679 
680  nameMass="histDstarMass_";
681  nameMass+=i+1;
682  nameSgn="histDstarSgn_";
683  nameSgn+=i+1;
684  nameBkg="histDstarBkg_";
685  nameBkg+=i+1;
686 
687  if (i==-2) {
688  nameMass="histDstarMass";
689  nameSgn="histDstarSgn";
690  nameBkg="histDstarBkg";
691  }
692 
693  TH1F* spectrumDstarMass = new TH1F(nameMass.Data(),"D^{*} invariant mass; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
694  TH1F* spectrumDstarSgn = new TH1F(nameSgn.Data(), "D^{*} Signal invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
695  TH1F* spectrumDstarBkg = new TH1F(nameBkg.Data(), "D^{*} Background invariant mass - MC; M(D^{*}) [GeV/c^{2}]; Entries",200,1.9,2.1);
696 
697  nameMass="histSideBandMass_";
698  nameMass+=i+1;
699  if (i==-2) {
700  nameMass="histSideBandMass";
701  }
702 
703  TH1F* spectrumSideBandMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} sideband mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
704 
705  nameMass="histWrongSignMass_";
706  nameMass+=i+1;
707  if (i==-2) {
708  nameMass="histWrongSignMass";
709  }
710 
711  TH1F* spectrumWrongSignMass = new TH1F(nameMass.Data(),"D^{*}-D^{0} wrongsign mass; M(D^{*}) [GeV/c^{2}]; Entries",200,0.1,0.2);
712 
713 
714  spectrumMass->Sumw2();
715  spectrumSgn->Sumw2();
716  spectrumBkg->Sumw2();
717 
718  spectrumMass->SetLineColor(6);
719  spectrumSgn->SetLineColor(2);
720  spectrumBkg->SetLineColor(4);
721 
722  spectrumMass->SetMarkerStyle(20);
723  spectrumSgn->SetMarkerStyle(20);
724  spectrumBkg->SetMarkerStyle(20);
725  spectrumMass->SetMarkerSize(0.6);
726  spectrumSgn->SetMarkerSize(0.6);
727  spectrumBkg->SetMarkerSize(0.6);
728  spectrumMass->SetMarkerColor(6);
729  spectrumSgn->SetMarkerColor(2);
730  spectrumBkg->SetMarkerColor(4);
731 
732  spectrumD0Mass->Sumw2();
733  spectrumD0Sgn->Sumw2();
734  spectrumD0Bkg->Sumw2();
735 
736  spectrumD0Mass->SetLineColor(6);
737  spectrumD0Sgn->SetLineColor(2);
738  spectrumD0Bkg->SetLineColor(4);
739 
740  spectrumD0Mass->SetMarkerStyle(20);
741  spectrumD0Sgn->SetMarkerStyle(20);
742  spectrumD0Bkg->SetMarkerStyle(20);
743  spectrumD0Mass->SetMarkerSize(0.6);
744  spectrumD0Sgn->SetMarkerSize(0.6);
745  spectrumD0Bkg->SetMarkerSize(0.6);
746  spectrumD0Mass->SetMarkerColor(6);
747  spectrumD0Sgn->SetMarkerColor(2);
748  spectrumD0Bkg->SetMarkerColor(4);
749 
750  spectrumDstarMass->Sumw2();
751  spectrumDstarSgn->Sumw2();
752  spectrumDstarBkg->Sumw2();
753 
754  spectrumDstarMass->SetLineColor(6);
755  spectrumDstarSgn->SetLineColor(2);
756  spectrumDstarBkg->SetLineColor(4);
757 
758  spectrumDstarMass->SetMarkerStyle(20);
759  spectrumDstarSgn->SetMarkerStyle(20);
760  spectrumDstarBkg->SetMarkerStyle(20);
761  spectrumDstarMass->SetMarkerSize(0.6);
762  spectrumDstarSgn->SetMarkerSize(0.6);
763  spectrumDstarBkg->SetMarkerSize(0.6);
764  spectrumDstarMass->SetMarkerColor(6);
765  spectrumDstarSgn->SetMarkerColor(2);
766  spectrumDstarBkg->SetMarkerColor(4);
767 
768  spectrumSideBandMass->Sumw2();
769  spectrumSideBandMass->SetLineColor(4);
770  spectrumSideBandMass->SetMarkerStyle(20);
771  spectrumSideBandMass->SetMarkerSize(0.6);
772  spectrumSideBandMass->SetMarkerColor(4);
773 
774  spectrumWrongSignMass->Sumw2();
775  spectrumWrongSignMass->SetLineColor(4);
776  spectrumWrongSignMass->SetMarkerStyle(20);
777  spectrumWrongSignMass->SetMarkerSize(0.6);
778  spectrumWrongSignMass->SetMarkerColor(4);
779 
780  TH1F* allMass = (TH1F*)spectrumMass->Clone();
781  TH1F* allSgn = (TH1F*)spectrumSgn->Clone();
782  TH1F* allBkg = (TH1F*)spectrumBkg->Clone();
783 
784  TH1F* pidMass = (TH1F*)spectrumMass->Clone();
785  TH1F* pidSgn = (TH1F*)spectrumSgn->Clone();
786  TH1F* pidBkg = (TH1F*)spectrumBkg->Clone();
787 
788  fOutputAll->Add(allMass);
789  fOutputAll->Add(allSgn);
790  fOutputAll->Add(allBkg);
791  fAllhist[i+2+((fNPtBins+2)*kDeltaMass)]=allMass;
792  fAllhist[i+2+((fNPtBins+2)*kDeltaSgn)]=allSgn;
793  fAllhist[i+2+((fNPtBins+2)*kDeltaBkg)]=allBkg;
794 
795  fOutputPID->Add(pidMass);
796  fOutputPID->Add(pidSgn);
797  fOutputPID->Add(pidBkg);
798  fPIDhist[i+2+((fNPtBins+2)*kDeltaMass)]=pidMass;
799  fPIDhist[i+2+((fNPtBins+2)*kDeltaSgn)]=pidSgn;
800  fPIDhist[i+2+((fNPtBins+2)*kDeltaBkg)]=pidBkg;
801 
802  TH1F* allD0Mass = (TH1F*)spectrumD0Mass->Clone();
803  TH1F* allD0Sgn = (TH1F*)spectrumD0Sgn->Clone();
804  TH1F* allD0Bkg = (TH1F*)spectrumD0Bkg->Clone();
805 
806  TH1F* pidD0Mass = (TH1F*)spectrumD0Mass->Clone();
807  TH1F* pidD0Sgn = (TH1F*)spectrumD0Sgn->Clone();
808  TH1F* pidD0Bkg = (TH1F*)spectrumD0Bkg->Clone();
809 
810  fOutputAll->Add(allD0Mass);
811  fOutputAll->Add(allD0Sgn);
812  fOutputAll->Add(allD0Bkg);
813  fAllhist[i+2+((fNPtBins+2)*kDzMass)]=allD0Mass;
814  fAllhist[i+2+((fNPtBins+2)*kDzSgn)]=allD0Sgn;
815  fAllhist[i+2+((fNPtBins+2)*kDzBkg)]=allD0Bkg;
816 
817  fOutputPID->Add(pidD0Mass);
818  fOutputPID->Add(pidD0Sgn);
819  fOutputPID->Add(pidD0Bkg);
820  fPIDhist[i+2+((fNPtBins+2)*kDzMass)]=pidD0Mass;
821  fPIDhist[i+2+((fNPtBins+2)*kDzSgn)]=pidD0Sgn;
822  fPIDhist[i+2+((fNPtBins+2)*kDzBkg)]=pidD0Bkg;
823 
824  TH1F* allDstarMass = (TH1F*)spectrumDstarMass->Clone();
825  TH1F* allDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
826  TH1F* allDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
827 
828  TH1F* pidDstarMass = (TH1F*)spectrumDstarMass->Clone();
829  TH1F* pidDstarSgn = (TH1F*)spectrumDstarSgn->Clone();
830  TH1F* pidDstarBkg = (TH1F*)spectrumDstarBkg->Clone();
831 
832  fOutputAll->Add(allDstarMass);
833  fOutputAll->Add(allDstarSgn);
834  fOutputAll->Add(allDstarBkg);
835  fAllhist[i+2+((fNPtBins+2)*kDstarMass)]=allDstarMass;
836  fAllhist[i+2+((fNPtBins+2)*kDstarSgn)]=allDstarSgn;
837  fAllhist[i+2+((fNPtBins+2)*kDstarBkg)]=allDstarBkg;
838 
839  fOutputPID->Add(pidDstarMass);
840  fOutputPID->Add(pidDstarSgn);
841  fOutputPID->Add(pidDstarBkg);
842  fPIDhist[i+2+((fNPtBins+2)*kDstarMass)]=pidDstarMass;
843  fPIDhist[i+2+((fNPtBins+2)*kDstarSgn)]=pidDstarSgn;
844  fPIDhist[i+2+((fNPtBins+2)*kDstarBkg)]=pidDstarBkg;
845 
846  TH1F* allSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
847  TH1F* pidSideBandMass = (TH1F*)spectrumSideBandMass->Clone();
848 
849  fOutputAll->Add(allSideBandMass);
850  fOutputPID->Add(pidSideBandMass);
851  fAllhist[i+2+((fNPtBins+2)*kSideBandMass)]=allSideBandMass;
852  fPIDhist[i+2+((fNPtBins+2)*kSideBandMass)]=pidSideBandMass;
853 
854  TH1F* allWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
855  TH1F* pidWrongSignMass = (TH1F*)spectrumWrongSignMass->Clone();
856 
857  fOutputAll->Add(allWrongSignMass);
858  fOutputPID->Add(pidWrongSignMass);
859  fAllhist[i+2+((fNPtBins+2)*kWrongSignMass)]=allWrongSignMass;
860  fPIDhist[i+2+((fNPtBins+2)*kWrongSignMass)]=pidWrongSignMass;
861 
862  }
863 
864  // pt spectra
865  nameMass="ptMass";
866  nameSgn="ptSgn";
867  nameBkg="ptBkg";
868 
869  TH1F* ptspectrumMass = new TH1F(nameMass.Data(),"D^{*} p_{T}; p_{T} [GeV]; Entries",400,0,50);
870  TH1F* ptspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal p_{T} - MC; p_{T} [GeV]; Entries",400,0,50);
871  TH1F* ptspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background p_{T} - MC; p_{T} [GeV]; Entries",400,0,50);
872 
873  ptspectrumMass->Sumw2();
874  ptspectrumSgn->Sumw2();
875  ptspectrumBkg->Sumw2();
876 
877  ptspectrumMass->SetLineColor(6);
878  ptspectrumSgn->SetLineColor(2);
879  ptspectrumBkg->SetLineColor(4);
880 
881  ptspectrumMass->SetMarkerStyle(20);
882  ptspectrumSgn->SetMarkerStyle(20);
883  ptspectrumBkg->SetMarkerStyle(20);
884  ptspectrumMass->SetMarkerSize(0.6);
885  ptspectrumSgn->SetMarkerSize(0.6);
886  ptspectrumBkg->SetMarkerSize(0.6);
887  ptspectrumMass->SetMarkerColor(6);
888  ptspectrumSgn->SetMarkerColor(2);
889  ptspectrumBkg->SetMarkerColor(4);
890 
891  TH1F* ptallMass = (TH1F*)ptspectrumMass->Clone();
892  TH1F* ptallSgn = (TH1F*)ptspectrumSgn->Clone();
893  TH1F* ptallBkg = (TH1F*)ptspectrumBkg->Clone();
894 
895  TH1F* ptpidMass = (TH1F*)ptspectrumMass->Clone();
896  TH1F* ptpidSgn = (TH1F*)ptspectrumSgn->Clone();
897  TH1F* ptpidBkg = (TH1F*)ptspectrumBkg->Clone();
898 
899  fOutputAll->Add(ptallMass);
900  fOutputAll->Add(ptallSgn);
901  fOutputAll->Add(ptallBkg);
902  fAllhist[((fNPtBins+2)*kptMass)]=ptallMass;
903  fAllhist[((fNPtBins+2)*kptSgn)]=ptallSgn;
904  fAllhist[((fNPtBins+2)*kptBkg)]=ptallBkg;
905 
906  fOutputPID->Add(ptpidMass);
907  fOutputPID->Add(ptpidSgn);
908  fOutputPID->Add(ptpidBkg);
909  fPIDhist[(fNPtBins+2)*kptMass]=ptpidMass;
910  fPIDhist[(fNPtBins+2)*kptSgn]=ptpidSgn;
911  fPIDhist[(fNPtBins+2)*kptBkg]=ptpidBkg;
912  // eta spectra
913  nameMass="etaMass";
914  nameSgn="etaSgn";
915  nameBkg="etaBkg";
916 
917  TH1F* etaspectrumMass = new TH1F(nameMass.Data(),"D^{*} #eta; #eta; Entries",200,-1,1);
918  TH1F* etaspectrumSgn = new TH1F(nameSgn.Data(), "D^{*} Signal #eta - MC; #eta; Entries",200,-1,1);
919  TH1F* etaspectrumBkg = new TH1F(nameBkg.Data(), "D^{*} Background #eta - MC; #eta; Entries",200,-1,1);
920 
921  etaspectrumMass->Sumw2();
922  etaspectrumSgn->Sumw2();
923  etaspectrumBkg->Sumw2();
924 
925  etaspectrumMass->SetLineColor(6);
926  etaspectrumSgn->SetLineColor(2);
927  etaspectrumBkg->SetLineColor(4);
928 
929  etaspectrumMass->SetMarkerStyle(20);
930  etaspectrumSgn->SetMarkerStyle(20);
931  etaspectrumBkg->SetMarkerStyle(20);
932  etaspectrumMass->SetMarkerSize(0.6);
933  etaspectrumSgn->SetMarkerSize(0.6);
934  etaspectrumBkg->SetMarkerSize(0.6);
935  etaspectrumMass->SetMarkerColor(6);
936  etaspectrumSgn->SetMarkerColor(2);
937  etaspectrumBkg->SetMarkerColor(4);
938 
939  TH1F* etaallMass = (TH1F*)etaspectrumMass->Clone();
940  TH1F* etaallSgn = (TH1F*)etaspectrumSgn->Clone();
941  TH1F* etaallBkg = (TH1F*)etaspectrumBkg->Clone();
942 
943  TH1F* etapidMass = (TH1F*)etaspectrumMass->Clone();
944  TH1F* etapidSgn = (TH1F*)etaspectrumSgn->Clone();
945  TH1F* etapidBkg = (TH1F*)etaspectrumBkg->Clone();
946 
947  fOutputAll->Add(etaallMass);
948  fOutputAll->Add(etaallSgn);
949  fOutputAll->Add(etaallBkg);
950  fAllhist[(fNPtBins+2)*ketaMass]=etaallMass;
951  fAllhist[(fNPtBins+2)*ketaSgn]=etaallSgn;
952  fAllhist[(fNPtBins+2)*ketaBkg]=etaallBkg;
953 
954  fOutputPID->Add(etapidMass);
955  fOutputPID->Add(etapidSgn);
956  fOutputPID->Add(etapidBkg);
957  fPIDhist[(fNPtBins+2)*ketaMass]=etapidMass;
958  fPIDhist[(fNPtBins+2)*ketaSgn]=etapidSgn;
959  fPIDhist[(fNPtBins+2)*ketaBkg]=etapidBkg;
960 
961  if (fDoDStarVsY){
962  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.);
963  fOutputPID->Add(deltamassVsyVsPtPID);
964  }
965  return;
966 }
967 //________________________________________________________________________
969  //
971  //
972 
973  if(!isSel) return;
974 
975  // D0 window
976  Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
977  Double_t invmassD0 = part->InvMassD0();
978 
979 
980  Int_t ptbin=cuts->PtBin(part->Pt());
981  Double_t pt = part->Pt();
982  Double_t eta = part->Eta();
983 
984  Double_t invmassDelta = part->DeltaInvMass();
985  Double_t invmassDstar = part->InvMassDstarKpipi();
986 
987  TString fillthis="";
988  Bool_t massInRange=kFALSE;
989 
990  Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
991 
992  // delta M(Kpipi)-M(Kpi)
993  if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<fPeakWindow) massInRange=kTRUE;
994 
995  if(fUseMCInfo) {
996  if(isDStar==1) {
997  histlist[ptbin+1+((fNPtBins+2)*kDzSgn)]->Fill(invmassD0);
998  histlist[(fNPtBins+2)*kDzSgn]->Fill(invmassD0);
999  histlist[ptbin+1+((fNPtBins+2)*kDstarSgn)]->Fill(invmassDstar);
1000  histlist[(fNPtBins+2)*kDstarSgn]->Fill(invmassDstar);
1001  histlist[ptbin+1+((fNPtBins+2)*kDeltaSgn)]->Fill(invmassDelta);
1002  histlist[(fNPtBins+2)*kDeltaSgn]->Fill(invmassDelta);
1003  if (massInRange) {
1004  histlist[(fNPtBins+2)*kptSgn]->Fill(pt);
1005  histlist[(fNPtBins+2)*ketaSgn]->Fill(eta);
1006  }
1007  }
1008  else {//background
1009  histlist[ptbin+1+((fNPtBins+2)*kDzBkg)]->Fill(invmassD0);
1010  histlist[(fNPtBins+2)*kDzBkg]->Fill(invmassD0);
1011  histlist[ptbin+1+((fNPtBins+2)*kDstarBkg)]->Fill(invmassDstar);
1012  histlist[(fNPtBins+2)*kDstarBkg]->Fill(invmassDstar);
1013  histlist[ptbin+1+((fNPtBins+2)*kDeltaBkg)]->Fill(invmassDelta);
1014  histlist[(fNPtBins+2)*kDeltaBkg]->Fill(invmassDelta);
1015  if (massInRange) {
1016  histlist[(fNPtBins+2)*kptBkg]->Fill(pt);
1017  histlist[(fNPtBins+2)*ketaBkg]->Fill(eta);
1018  }
1019  }
1020  }
1021  //no MC info, just cut selection
1022  histlist[ptbin+1+((fNPtBins+2)*kDzMass)]->Fill(invmassD0);
1023  histlist[(fNPtBins+2)*kDzMass]->Fill(invmassD0);
1024  histlist[ptbin+1+((fNPtBins+2)*kDstarMass)]->Fill(invmassDstar);
1025  histlist[(fNPtBins+2)*kDstarMass]->Fill(invmassDstar);
1026  histlist[ptbin+1+((fNPtBins+2)*kDeltaMass)]->Fill(invmassDelta);
1027  histlist[(fNPtBins+2)*kDeltaMass]->Fill(invmassDelta);
1028 
1029  if (massInRange) {
1030  histlist[(fNPtBins+2)*kptMass]->Fill(pt);
1031  histlist[(fNPtBins+2)*ketaMass]->Fill(eta);
1032  }
1033 
1034  return;
1035 }
1036 //______________________________ side band background for D*___________________________________
1038 
1042 
1043  if(!isSel) return;
1044 
1045  Int_t ptbin=cuts->PtBin(part->Pt());
1046 
1047  // select the side bands intervall
1048  Double_t invmassD0 = part->InvMassD0();
1049  if(TMath::Abs(invmassD0-1.865)>4*fD0Window && TMath::Abs(invmassD0-1.865)<8*fD0Window){
1050 
1051  // for pt and eta
1052  Double_t invmassDelta = part->DeltaInvMass();
1053 
1054  histlist[ptbin+1+((fNPtBins+2)*kSideBandMass)]->Fill(invmassDelta);
1055  histlist[(fNPtBins+2)*kSideBandMass]->Fill(invmassDelta);
1056 
1057 
1058  }
1059 }
1060 //________________________________________________________________________________________________________________
1062  //
1064  //
1065  Int_t ptbin=cuts->PtBin(part->Pt());
1066 
1067  Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1068  Double_t invmassD0 = part->InvMassD0();
1069  if (TMath::Abs(invmassD0-mPDGD0)>fD0Window) return;
1070 
1071  AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)part->Get2Prong();
1072 
1073  Int_t okDzWrongSign;
1074  Double_t wrongMassD0=0.;
1075 
1076  Int_t isSelected=cuts->IsSelected(part,AliRDHFCuts::kCandidate); //selected
1077  if (!isSelected){
1078  return;
1079  }
1080 
1081  okDzWrongSign = 1;
1082 
1083  //if is D*+ than assume D0bar
1084  if(part->Charge()>0 && (isSelected ==1)) {
1085  okDzWrongSign = 0;
1086  }
1087 
1088  // assign the wrong mass in case the cuts return both D0 and D0bar
1089  if(part->Charge()>0 && (isSelected ==3)) {
1090  okDzWrongSign = 0;
1091  }
1092 
1093  //wrong D0 inv mass
1094  if(okDzWrongSign!=0){
1095  wrongMassD0 = theD0particle->InvMassD0();
1096  }else if(okDzWrongSign==0){
1097  wrongMassD0 = theD0particle->InvMassD0bar();
1098  }
1099 
1100  if(TMath::Abs(wrongMassD0-1.865)<fD0Window){
1101 
1102  // wrong D* inv mass
1103  Double_t e[3];
1104  if (part->Charge()>0){
1105  e[0]=theD0particle->EProng(0,321);
1106  e[1]=theD0particle->EProng(1,211);
1107  }else{
1108  e[0]=theD0particle->EProng(0,211);
1109  e[1]=theD0particle->EProng(1,321);
1110  }
1111  e[2]=part->EProng(0,211);
1112 
1113  Double_t esum = e[0]+e[1]+e[2];
1114  Double_t pds = part->P();
1115 
1116  Double_t wrongMassDstar = TMath::Sqrt(esum*esum-pds*pds);
1117 
1118  TString fillthis="";
1119  fillthis="histWrongSignMass_";
1120  fillthis+=ptbin;
1121  ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);
1122  fillthis="histWrongSignMass";
1123  ((TH1F*)(listout->FindObject(fillthis)))->Fill(wrongMassDstar-wrongMassD0);
1124 
1125  }
1126 }
1127 
1128 //-------------------------------------------------------------------------------
1129 Int_t AliAnalysisTaskSEDStarSpectra::CheckOrigin(TClonesArray* arrayMC, const AliAODMCParticle *mcPartCandidate) const {
1130  //
1131  // checking whether the mother of the particles come from a charm or a bottom quark
1132  //
1133 
1134  Int_t pdgGranma = 0;
1135  Int_t mother = 0;
1136  mother = mcPartCandidate->GetMother();
1137  Int_t istep = 0;
1138  Int_t abspdgGranma =0;
1139  Bool_t isFromB=kFALSE;
1140  while (mother >0 ){
1141  istep++;
1142  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
1143  if (mcGranma){
1144  pdgGranma = mcGranma->GetPdgCode();
1145  abspdgGranma = TMath::Abs(pdgGranma);
1146  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
1147  isFromB=kTRUE;
1148  }
1149  mother = mcGranma->GetMother();
1150  }else{
1151  AliError("Failed casting the mother particle!");
1152  break;
1153  }
1154  }
1155 
1156  if(isFromB) return 5;
1157  else return 4;
1158 }
1159 //-------------------------------------------------------------------------------------
1160 Float_t AliAnalysisTaskSEDStarSpectra::GetTrueImpactParameterD0(const AliAODMCHeader *mcHeader, TClonesArray* arrayMC, const AliAODMCParticle *partDp) const {
1162 
1163  Double_t vtxTrue[3];
1164  mcHeader->GetVertex(vtxTrue);
1165  Double_t origD[3];
1166  partDp->XvYvZv(origD);
1167  Short_t charge=partDp->Charge();
1168  Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
1169  Int_t labelFirstDau = partDp->GetDaughter(0);
1170 
1171  Int_t nDau=partDp->GetNDaughters();
1172 
1173  Int_t theDau=0;
1174  if(nDau==2){
1175  for(Int_t iDau=0; iDau<2; iDau++){
1176  Int_t ind = labelFirstDau+iDau;
1177  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
1178  if(!part){
1179  AliError("Daughter particle not found in MC array");
1180  return 99999.;
1181  }
1182  Int_t pdgCode=TMath::Abs(part->GetPdgCode());
1183  if(pdgCode==211 || pdgCode==321){
1184  pXdauTrue[theDau]=part->Px();
1185  pYdauTrue[theDau]=part->Py();
1186  pZdauTrue[theDau]=part->Pz();
1187  ++theDau;
1188  }
1189  }
1190  }
1191  if(theDau!=2){
1192  AliError("Wrong number of decay prongs");
1193  return 99999.;
1194  }
1195 
1196  Double_t d0dummy[3]={0.,0.,0.};
1197  AliAODRecoDecayHF aodD0MC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
1198  return aodD0MC.ImpParXY();
1199 
1200 }
1201 //______________________________________________________-
1204 
1205  Int_t nbins[3]={400,200,fNImpParBins};
1206  Double_t xmin[3]={1.75,0.,fLowerImpPar};
1207  Double_t xmax[3]={1.98,20.,fHigherImpPar};
1208 
1209  fHistMassPtImpParTCDs[0]=new THnSparseF("hMassPtImpParAll",
1210  "Mass vs. pt vs.imppar - All",
1211  3,nbins,xmin,xmax);
1212  fHistMassPtImpParTCDs[1]=new THnSparseF("hMassPtImpParPrompt",
1213  "Mass vs. pt vs.imppar - promptD",
1214  3,nbins,xmin,xmax);
1215  fHistMassPtImpParTCDs[2]=new THnSparseF("hMassPtImpParBfeed",
1216  "Mass vs. pt vs.imppar - DfromB",
1217  3,nbins,xmin,xmax);
1218  fHistMassPtImpParTCDs[3]=new THnSparseF("hMassPtImpParTrueBfeed",
1219  "Mass vs. pt vs.true imppar -DfromB",
1220  3,nbins,xmin,xmax);
1221  fHistMassPtImpParTCDs[4]=new THnSparseF("hMassPtImpParBkg",
1222  "Mass vs. pt vs.imppar - backgr.",
1223  3,nbins,xmin,xmax);
1224 
1225  for(Int_t i=0; i<5;i++){
1226  fOutput->Add(fHistMassPtImpParTCDs[i]);
1227  }
1228 }
1229 
Int_t charge
AliAODTrack * GetProng(AliVEvent *event, AliAODRecoDecayHF *rd, Int_t iprong)
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)
virtual Int_t PreSelect(TObjArray aodtracks)
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:317
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)
UShort_t GetProngID(Int_t ip) const
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:210
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:250
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