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