AliPhysics  d565ceb (d565ceb)
AliAnalysisTaskSESignificance.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2008, 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  * appear 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 
19 //
20 // AliAnalysisTaskSESignificane to calculate effects on
21 // significance of D mesons cut
22 // Authors: G. Ortona, ortona@to.infn.it
23 // F. Prino, prino@to.infn.it
24 // Renu Bala, bala@to.infn.it
25 // Chiara Bianchin, cbianchi@pd.infn.it
27 
28 #include <Riostream.h>
29 #include <TClonesArray.h>
30 #include <TCanvas.h>
31 #include <TList.h>
32 #include <TFile.h>
33 #include <TString.h>
34 #include <TH1F.h>
35 #include <TDatabasePDG.h>
36 
37 #include <AliLog.h>
38 #include "AliAnalysisManager.h"
39 #include "AliAODHandler.h"
40 #include "AliAODEvent.h"
41 #include "AliAODVertex.h"
42 #include "AliAODTrack.h"
43 #include "AliAODMCHeader.h"
44 #include "AliAODMCParticle.h"
46 #include "AliAODRecoDecayHF.h"
49 #include "AliAODRecoCascadeHF.h"
50 
51 #include "AliAnalysisTaskSE.h"
53 #include "AliRDHFCutsD0toKpipipi.h"
54 #include "AliRDHFCutsDstoKKpi.h"
56 #include "AliRDHFCutsD0toKpi.h"
57 #include "AliRDHFCutsLctopKpi.h"
58 #include "AliMultiDimVector.h"
59 
61 
62 using std::cout;
63 using std::endl;
64 
68 
69 //________________________________________________________________________
72  fOutput(0),
73  fCutList(0),
74  fHistNEvents(0),
75  fUpmasslimit(1.965),
76  fLowmasslimit(1.765),
77  fRDCuts(0),
78  fNPtBins(0),
79  fAODProtection(1),
80  fReadMC(kFALSE),
81  fUseSelBit(kFALSE),
82  fBFeedDown(kBoth),
83  fDecChannel(0),
84  fPDGmother(0),
85  fNProngs(0),
86  fBranchName(""),
87  fSelectionlevel(0),
88  fNVars(0),
89  fNBins(100),
90  fPartOrAndAntiPart(0),
91  fDsChannel(0)
92 {
93  // Default constructor
94  SetPDGCodes();
96 
97  for(Int_t i=0;i<4;i++) fPDGdaughters[i]=0;
98  for(Int_t i=0;i<2;i++) {fPDGD0ToKpi[i]=0;fPDGDStarToD0pi[i]=0;}
99  for(Int_t i=0;i<kMaxCutVar;i++) fVars[i]=0.;
100  for(Int_t i=0;i<kMaxNHist;i++) {
101  fMassHist[i]=0;
102  fSigHist[i]=0;
103  fBkgHist[i]=0;
104  fRflHist[i]=0;
105  }
106 
107 }
108 
109 //________________________________________________________________________
110 AliAnalysisTaskSESignificance::AliAnalysisTaskSESignificance(const char *name, TList* listMDV,AliRDHFCuts *rdCuts,Int_t decaychannel,Int_t selectionlevel):
111  AliAnalysisTaskSE(name),
112  fOutput(0),
113  fCutList(listMDV),
114  fHistNEvents(0),
115  fUpmasslimit(0),
116  fLowmasslimit(0),
117  fRDCuts(rdCuts),
118  fNPtBins(0),
119  fAODProtection(1),
120  fReadMC(kFALSE),
121  fUseSelBit(kFALSE),
122  fBFeedDown(kBoth),
123  fDecChannel(decaychannel),
124  fPDGmother(0),
125  fNProngs(0),
126  fBranchName(""),
127  fSelectionlevel(selectionlevel),
128  fNVars(0),
129  fNBins(100),
131  fDsChannel(0)
132 {
133 
134  SetPDGCodes();
136  if (fDecChannel!=2) SetMassLimits(0.15,fPDGmother); //check range
137  else {
138  Float_t min = 0.13;
139  Float_t max = 0.19;
140  SetMassLimits(min, max);
141  }
143 
145  if(fNVars>kMaxCutVar) AliFatal(Form("Too large number of cut variables, maximum is %d",kMaxCutVar));
146 
147  if(fDebug>1)fRDCuts->PrintAll();
148  // Output slot #1 writes into a TList container
149  DefineOutput(1,TList::Class()); //My private output
150  DefineOutput(2,TList::Class());
151  DefineOutput(3,AliRDHFCuts::Class()); //class of the cuts
153 }
154 
155  //________________________________________________________________________
157 {
158  // Destructor
159  if (fOutput) {
160  delete fOutput;
161  fOutput = 0;
162  }
163  if (fCutList) {
164  delete fCutList;
165  fCutList = 0;
166  }
167  if(fHistNEvents){
168  delete fHistNEvents;
169  fHistNEvents=0;
170  }
171 /*
172  if(fRDCuts) {
173  delete fRDCuts;
174  fRDCuts= 0;
175  }
176 */
177 
178 }
179 //_________________________________________________________________
181  // sets channel dependent quantities
182 
183  switch(fDecChannel){
184  case 0:
185  //D+
186  fPDGmother=411;
187  fNProngs=3;
188  fPDGdaughters[0]=211;//pi
189  fPDGdaughters[1]=321;//K
190  fPDGdaughters[2]=211;//pi
191  fPDGdaughters[3]=0; //empty
192  fBranchName="Charm3Prong";
193  break;
194  case 1:
195  //D0
196  fPDGmother=421;
197  fNProngs=2;
198  fPDGdaughters[0]=211;//pi
199  fPDGdaughters[1]=321;//K
200  fPDGdaughters[2]=0; //empty
201  fPDGdaughters[3]=0; //empty
202  fBranchName="D0toKpi";
203  break;
204  case 2:
205  //D*
206  fPDGmother=413;
207  fNProngs=3;
208  fPDGdaughters[1]=211;//pi
209  fPDGdaughters[0]=321;//K
210  fPDGdaughters[2]=211;//pi (soft?)
211  fPDGdaughters[3]=0; //empty
212  fBranchName="Dstar";
213  break;
214  case 3:
215  //Ds
216  fPDGmother=431;
217  fNProngs=3;
218  fPDGdaughters[0]=321;//K
219  fPDGdaughters[1]=321;//K
220  fPDGdaughters[2]=211;//pi
221  fPDGdaughters[3]=0; //empty
222  fBranchName="Charm3Prong";
223  break;
224  case 4:
225  //D0 in 4 prongs
226  fPDGmother=421;
227  fNProngs=4;
228  fPDGdaughters[0]=321;
229  fPDGdaughters[1]=211;
230  fPDGdaughters[2]=211;
231  fPDGdaughters[3]=211;
232  fBranchName="Charm4Prong";
233  break;
234  case 5:
235  //Lambda_c
236  fPDGmother=4122;
237  fNProngs=3;
238  fPDGdaughters[0]=2212;//p
239  fPDGdaughters[1]=321;//K
240  fPDGdaughters[2]=211;//pi
241  fPDGdaughters[3]=0; //empty
242  fBranchName="Charm3Prong";
243  break;
244  }
245 }
246 
247 //_________________________________________________________________
249 
250  Bool_t result = kTRUE;
251 
252  const Int_t nvars=fRDCuts->GetNVars();//ForOpt();
253  //Float_t *vars = new Float_t[nvars];
254  Bool_t *varsforopt = fRDCuts->GetVarsForOpt();
255  Bool_t *uppervars = fRDCuts->GetIsUpperCut();
256  TString *names = fRDCuts->GetVarNames();
257 
258  for(Int_t i=0;i<fNPtBins;i++){
259  TString mdvname=Form("multiDimVectorPtBin%d",i);
260  Int_t ic=0;
261 
262  for(Int_t ivar=0;ivar<nvars;ivar++){
263  if(varsforopt[ivar]){
264  Float_t min = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetMinLimit(ic);
265  Float_t max = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetMaxLimit(ic);
266  if(min==max){
267  AliFatal(Form("tight and loose cut for optimization variable number %d are the same in ptbin %d\n",ic,i));
268  return kFALSE;
269  }
270  Bool_t lowermdv = ((AliMultiDimVector*)fCutList->FindObject(mdvname.Data()))->GetGreaterThan(ic);
271  if(uppervars[ivar]&&lowermdv){
272  AliFatal(Form("%s is declared as uppercut, but as been given tighter cut larger then loose cut in ptbin %d \n ---please check your cuts \n ",names[ivar].Data(),i));
273  return kFALSE;
274  }
275  if(!uppervars[ivar]&&!lowermdv){
276  AliFatal(Form("%s is declared as lower cut, but as been given tighter cut smaller then loose cut in ptbin %d \n ---please check your cuts \n",names[ivar].Data(),i));
277  return kFALSE;
278  }
279  ic++;
280  }
281  }
282  }
283  return result;
284 }
285 //_________________________________________________________________
287  if(fReadMC)fBFeedDown=flagB;
288  else {
289  AliInfo("B feed down not allowed without MC info\n");
291  }
292 }
293 //_________________________________________________________________
295  Float_t mass=0;
296  Int_t abspdg=TMath::Abs(pdg);
297  mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
298  fUpmasslimit = mass+range;
299  fLowmasslimit = mass-range;
300 }
301 //_________________________________________________________________
303  if(uplimit>lowlimit)
304  {
305  fUpmasslimit = uplimit;
306  fLowmasslimit = lowlimit;
307  }
308 }
309 
310 
311 
312 //________________________________________________________________________
314 {
315  // Initialization
316 
317  if(fDebug > 1) printf("AnalysisTaskSESignificance::Init() \n");
318 
319  switch(fDecChannel){
320  case 0:
321  {
322  AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCuts)));
323  // Post the data
324  PostData(3,copycut);
325  }
326  break;
327  case 1:
328  {
329  AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCuts)));
330  // Post the data
331  PostData(3,copycut);
332  }
333  break;
334  case 2:
335  {
336  AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCuts)));
337  // Post the data
338  PostData(3,copycut);
339  }
340  break;
341  case 3:
342  {
343  AliRDHFCutsDstoKKpi* copycut=new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fRDCuts)));
344  // Post the data
345  PostData(3,copycut);
346  }
347  break;
348  case 4:
349  {
350  AliRDHFCutsD0toKpipipi* copycut=new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fRDCuts)));
351  // Post the data
352  PostData(3,copycut);
353  }
354  break;
355  case 5:
356  {
357  AliRDHFCutsLctopKpi* copycut=new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fRDCuts)));
358  // Post the data
359  PostData(3,copycut);
360  }
361  break;
362 
363  default:
364  return;
365  }
366 
367  TList *mdvList = new TList();
368  mdvList->SetOwner();
369  mdvList = fCutList;
370 
371  PostData(2,mdvList);
372 
373 
374 }
375 //________________________________________________________________________
377 {
378  // Create the output container
379 
380  if(fDebug > 1) printf("AnalysisTaskSESignificance::UserCreateOutputObjects() \n");
381 
382  // Several histograms are more conveniently managed in a TList
383  fOutput = new TList();
384  fOutput->SetOwner();
385  fOutput->SetName("OutputHistos");
386 
387  //same number of steps in each multiDimVectorPtBin%d !
388  Int_t nHist=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
389  cout<<"ncells = "<<nHist<<" n ptbins = "<<fNPtBins<<endl;
390  nHist=nHist*fNPtBins;
391  cout<<"Total = "<<nHist<<endl;
392  for(Int_t i=0;i<nHist;i++){
393 
394  TString hisname;
395  TString signame;
396  TString bkgname;
397  TString rflname;
398  TString title;
399 
400  hisname.Form("hMass_%d",i);
401  signame.Form("hSig_%d",i);
402  bkgname.Form("hBkg_%d",i);
403  rflname.Form("hRfl_%d",i);
404 
405  title.Form("Invariant mass;M[GeV/c^{2}];Entries");
406 
407  fMassHist[i]=new TH1F(hisname.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
408  fMassHist[i]->Sumw2();
409  fOutput->Add(fMassHist[i]);
410 
411  if(fReadMC){
412  fSigHist[i]=new TH1F(signame.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
413  fSigHist[i]->Sumw2();
414  fOutput->Add(fSigHist[i]);
415 
416  fBkgHist[i]=new TH1F(bkgname.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
417  fBkgHist[i]->Sumw2();
418  fOutput->Add(fBkgHist[i]);
419 
421  fRflHist[i]=new TH1F(rflname.Data(),title.Data(),fNBins,fLowmasslimit,fUpmasslimit);
422  fRflHist[i]->Sumw2();
423  fOutput->Add(fRflHist[i]);
424  }
425  }
426  }
427 
428  fHistNEvents=new TH1F("fHistNEvents","Number of AODs scanned",8,-0.5,7.5);
429  fHistNEvents->GetXaxis()->SetBinLabel(1,"nEventsAnal");
430  fHistNEvents->GetXaxis()->SetBinLabel(2,"nEvSelected (vtx)");
431  fHistNEvents->GetXaxis()->SetBinLabel(3,"nCandidatesSelected");
432  fHistNEvents->GetXaxis()->SetBinLabel(4,"nTotEntries Mass hists");
433  fHistNEvents->GetXaxis()->SetBinLabel(5,"Pile-up Rej");
434  fHistNEvents->GetXaxis()->SetBinLabel(6,"N. of 0SMH");
435  if(fReadMC){
436  fHistNEvents->GetXaxis()->SetBinLabel(7,"MC Cand from c");
437  fHistNEvents->GetXaxis()->SetBinLabel(8,"MC Cand from b");
438  } else{
439  fHistNEvents->GetXaxis()->SetBinLabel(7,"N candidates");
440  }
441  fHistNEvents->GetXaxis()->SetNdivisions(1,kFALSE);
442  fOutput->Add(fHistNEvents);
443 
444 
445  PostData(1,fOutput);
446  return;
447 }
448 
449 //________________________________________________________________________
451 {
452  // Execute analysis for current event:
453  // heavy flavor candidates association to MC truth
454 
455  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
456  if(fAODProtection>=0){
457  // Protection against different number of events in the AOD and deltaAOD
458  // In case of discrepancy the event is rejected.
459  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
460  if (matchingAODdeltaAODlevel<0 || (matchingAODdeltaAODlevel==0 && fAODProtection==1)) {
461  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
462  return;
463  }
464  }
465  if(fDebug>2) printf("Analysing decay %d\n",fDecChannel);
466  // Post the data already here
467  PostData(1,fOutput);
468  TClonesArray *arrayProng =0;
469 
470  if(!aod && AODEvent() && IsStandardAOD()) {
471  // In case there is an AOD handler writing a standard AOD, use the AOD
472  // event in memory rather than the input (ESD) event.
473  aod = dynamic_cast<AliAODEvent*> (AODEvent());
474  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
475  // have to taken from the AOD event hold by the AliAODExtension
476  AliAODHandler* aodHandler = (AliAODHandler*)
477  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
478  if(aodHandler->GetExtensions()) {
479 
480  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
481  AliAODEvent *aodFromExt = ext->GetAOD();
482  arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
483 
484  }
485  } else if(aod) {
486  arrayProng=(TClonesArray*)aod->GetList()->FindObject(fBranchName.Data());
487  }
488  if(!aod || !arrayProng) {
489  AliError("AliAnalysisTaskSESignificance::UserExec:Branch not found!\n");
490  return;
491  }
492 
493  // fix for temporary bug in ESDfilter
494  // the AODs with null vertex pointer didn't pass the PhysSel
495  if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
496  TClonesArray *arrayMC=0;
497  AliAODMCHeader *mcHeader=0;
498 
499  // load MC particles
500  if(fReadMC){
501 
502  arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
503  if(!arrayMC) {
504  AliError("AliAnalysisTaskSESignificance::UserExec:MC particles branch not found!\n");
505  return;
506  }
507 
508  // load MC header
509  mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
510  if(!mcHeader) {
511  AliError("AliAnalysisTaskSESignificance::UserExec:MC header branch not found!\n");
512  return;
513  }
514  }
515 
516  fHistNEvents->Fill(0); // count event
517 
518  AliAODRecoDecayHF *d=0;
519 
520 
521 
522  Int_t nHistpermv=((AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0"))->GetNTotCells();
523  Int_t nProng = arrayProng->GetEntriesFast();
524  if(fDebug>1) printf("Number of D2H: %d\n",nProng);
525 
526  // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
527  TString trigclass=aod->GetFiredTriggerClasses();
528  if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fHistNEvents->Fill(5);
529 
530  if(fRDCuts->IsEventSelected(aod)) {
531  fHistNEvents->Fill(1);
532  }else{
533  if(fRDCuts->GetWhyRejection()==1) // rejected for pileup
534  fHistNEvents->Fill(4);
535  return;
536  }
537 
539 
540  for (Int_t iProng = 0; iProng < nProng; iProng++) {
541  fHistNEvents->Fill(6);
542  d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iProng);
543 
544  AliAODRecoCascadeHF* DStarToD0pi = NULL;
545  AliAODRecoDecayHF2Prong* D0Particle = NULL;
546  fPDGDStarToD0pi[0] = 421; fPDGDStarToD0pi[1] = 211;
547  fPDGD0ToKpi[0] = 321; fPDGD0ToKpi[1] = 211;
548 
549  Bool_t isSelBit=kTRUE;
550  if(fUseSelBit){
551  if(fDecChannel==0) {
553  }else{
554  if(fDecChannel==1) {
556  }else{
557  if(fDecChannel==2) {
559  }else{
560  if(fDecChannel==3) {
562  }else{
564  }
565  }
566  }
567  }
568  }
569  if(!isSelBit) continue;
570 
571  if(fDecChannel==0 || fDecChannel==3) {
572  if(!vHF->FillRecoCand(aod,(AliAODRecoDecayHF3Prong*)d))continue;
573  }
574  else if(fDecChannel == 1) {
575  if(!vHF->FillRecoCand(aod,(AliAODRecoDecayHF2Prong*)d))continue;
576  }
577  else if(fDecChannel == 2) {
578  if(!vHF->FillRecoCasc(aod,((AliAODRecoCascadeHF*)d),kTRUE))continue;
579  }
580 
581  if (fDecChannel==2) {
582  DStarToD0pi = (AliAODRecoCascadeHF*)arrayProng->At(iProng);
583  if (!DStarToD0pi->GetSecondaryVtx()) continue;
584  D0Particle = (AliAODRecoDecayHF2Prong*)DStarToD0pi->Get2Prong();
585  if (!D0Particle) continue;
586  }
587 
588  Bool_t isFidAcc = fRDCuts->IsInFiducialAcceptance(d->Pt(),d->Y(fPDGmother));
589  Int_t isSelected=fRDCuts->IsSelected(d,fSelectionlevel,aod);
590 
591  if(fReadMC && fBFeedDown!=kBoth && isSelected){
592  Int_t labD = d->MatchToMC(fPDGmother,arrayMC,fNProngs,fPDGdaughters);
593  if(labD>=0){
594  AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
595  Int_t pdgGranma = CheckOrigin(partD, arrayMC);
596  Int_t abspdgGranma = TMath::Abs(pdgGranma);
597  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
598  //feed down particle
599  AliDebug(2,Form("Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
600  fHistNEvents->Fill(7);
601  if(fBFeedDown==kCharmOnly) isSelected=kFALSE; //from beauty
602  }
603  else {
604  //prompt particle
605  fHistNEvents->Fill(6);
606  if(fBFeedDown==kBeautyOnly)isSelected=kFALSE;
607  }
608  }
609  }
610 
611  if(isSelected&&isFidAcc) {
612  fHistNEvents->Fill(2); // count selected with loosest cuts
613  if(fDebug>1) printf("+++++++Is Selected\n");
614 
615  Int_t nVals=0;
618  Int_t ptbin=fRDCuts->PtBin(d->Pt());
619  if(ptbin==-1) continue;
620  TString mdvname=Form("multiDimVectorPtBin%d",ptbin);
621  AliMultiDimVector* muvec=(AliMultiDimVector*)fCutList->FindObject(mdvname.Data());
622 
623  ULong64_t *addresses = muvec->GetGlobalAddressesAboveCuts(fVars,(Float_t)d->Pt(),nVals);
624  if(fDebug>1)printf("nvals = %d\n",nVals);
625  for(Int_t ivals=0;ivals<nVals;ivals++){
626  if(addresses[ivals]>=muvec->GetNTotCells()){
627  if (fDebug>1) printf("Overflow!!\n");
628  delete [] addresses;
629  return;
630  }
631 
632  fHistNEvents->Fill(3);
633 
634  //fill the histograms with the appropriate method
635  switch (fDecChannel){
636  case 0:
637  FillDplus(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
638  break;
639  case 1:
640  FillD02p(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
641  break;
642  case 2:
643  FillDstar(DStarToD0pi,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
644  break;
645  case 3:
646  if(isSelected&1){
647  FillDs(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected,1);
648  }
649  break;
650  case 4:
651  FillD04p(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
652  break;
653  case 5:
654  FillLambdac(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
655  break;
656  default:
657  break;
658  }
659 
660  }
661 
662  if (fDecChannel==3 && isSelected&2){
664  nVals=0;
666  delete [] addresses;
667  addresses = muvec->GetGlobalAddressesAboveCuts(fVars,(Float_t)d->Pt(),nVals);
668  if(fDebug>1)printf("nvals = %d\n",nVals);
669  for(Int_t ivals=0;ivals<nVals;ivals++){
670  if(addresses[ivals]>=muvec->GetNTotCells()){
671  if (fDebug>1) printf("Overflow!!\n");
672  delete [] addresses;
673  return;
674  }
675  FillDs(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected,0);
676 
677  }
678 
679  }
680 
681  delete [] addresses;
682  }// end if selected
683 
684  }
685 
686 
687  PostData(1,fOutput);
688  return;
689 }
690 
691 //***************************************************************************
692 
693 // Methods used in the UserExec
694 
695 
696 //********************************************************************************************
697 
698 //Methods to fill the istograms with MC information, one for each candidate
699 //NB: the implementation for each candidate is responsibility of the corresponding developer
700 
701 void AliAnalysisTaskSESignificance::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel){
702  //D+ channel
703  if(!isSel){
704  AliError("Candidate not selected\n");
705  return;
706  }
707 
708  if(fPartOrAndAntiPart*d->GetCharge()<0)return;
709  Int_t pdgdaughters[3] = {211,321,211};
710  Double_t mass=d->InvMass(3,(UInt_t*)pdgdaughters);
711 
712  fMassHist[index]->Fill(mass);
713 
714 
715  if(fReadMC){
716  Int_t lab=-1;
717  lab = d->MatchToMC(411,arrayMC,3,pdgdaughters);
718  if(lab>=0){ //signal
719  fSigHist[index]->Fill(mass);
720  } else{ //background
721  fBkgHist[index]->Fill(mass);
722  }
723  }
724 }
725 
726 
727 void AliAnalysisTaskSESignificance::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel){
728 
729  //D0->Kpi channel
730 
731  //mass histograms
732  Int_t pdgdaughtersD0[2]={211,321};//pi,K
733  Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
734  Double_t masses[2];
735  masses[0]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
736  masses[1]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
737 
738  if((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0) fMassHist[index]->Fill(masses[0]);
739  if(isSel>=2 && fPartOrAndAntiPart<=0) fMassHist[index]->Fill(masses[1]);
740 
741 
742 
743  //MC histograms
744  if(fReadMC){
745 
746  Int_t matchtoMC=-1;
747 
748  //D0
749  Int_t pdgdaughters[2];
750  pdgdaughters[0]=211;//pi
751  pdgdaughters[1]=321;//K
752 
753  matchtoMC = d->MatchToMC(fPDGmother,arrayMC,fNProngs,pdgdaughters);
754 
755  Int_t prongPdgPlus=421,prongPdgMinus=(-1)*421;
756  if((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0){ //D0
757  if(matchtoMC>=0){
758  AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
759  Int_t pdgMC = dMC->GetPdgCode();
760 
761  if(pdgMC==prongPdgPlus) fSigHist[index]->Fill(masses[0]);
762  else fRflHist[index]->Fill(masses[0]);
763 
764  } else fBkgHist[index]->Fill(masses[0]);
765 
766  }
767  if(isSel>=2 && fPartOrAndAntiPart<=0){ //D0bar
768  if(matchtoMC>=0){
769  AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
770  Int_t pdgMC = dMC->GetPdgCode();
771 
772  if(pdgMC==prongPdgMinus) fSigHist[index]->Fill(masses[1]);
773  else fRflHist[index]->Fill(masses[1]);
774  } else fBkgHist[index]->Fill(masses[1]);
775  }
776  }
777 }
778 
779 void AliAnalysisTaskSESignificance::FillDstar(AliAODRecoCascadeHF* dstarD0pi,TClonesArray *arrayMC,Int_t index,Int_t isSel){
780  //D* channel
781 
782  AliInfo("Dstar selected\n");
783 
784  Double_t mass = dstarD0pi->DeltaInvMass();
785 
786  if((isSel>0) && TMath::Abs(fPartOrAndAntiPart)>=0) fMassHist[index]->Fill(mass);
787 
788  if(fReadMC) {
789  Int_t matchtoMC = -1;
790  matchtoMC = dstarD0pi->MatchToMC(413,421,fPDGDStarToD0pi, fPDGD0ToKpi, arrayMC);
791 
792  Int_t prongPdgDStarPlus=413;//,prongPdgDStarMinus=(-1)*413;
793 
794  if ((isSel>1) && TMath::Abs(fPartOrAndAntiPart)>=0) {
795  //D*+
796  if(matchtoMC>=0) {
797  AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
798  Int_t pdgMC = dMC->GetPdgCode();
799 
800  if (pdgMC==prongPdgDStarPlus) fSigHist[index]->Fill(mass);
801  else {
802  dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
803  mass = dstarD0pi->DeltaInvMass();
804  fRflHist[index]->Fill(mass);
805  dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
806  }
807  }
808  else fBkgHist[index]->Fill(mass);
809  }
810  }
811 }
812 
813 
814 void AliAnalysisTaskSESignificance::FillDs(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel,Int_t optDecay){
815 
816  // Fill Ds histos
817 
818 
819  Int_t pdgDsKKpi[3]={321,321,211};//K,K,pi
820  Int_t pdgDspiKK[3]={211,321,321};//pi,K,K
821  Double_t masses[2];
822  masses[0]=d->InvMass(fNProngs,(UInt_t*)pdgDsKKpi); //Ds
823  masses[1]=d->InvMass(fNProngs,(UInt_t*)pdgDspiKK); //Dsbar
824 
825  Int_t labDs=-1;
826  if(fReadMC){
827  labDs = d->MatchToMC(431,arrayMC,3,pdgDsKKpi);
828  }
829 
830  Int_t isKKpi=isSel&1;
831  Int_t ispiKK=isSel&2;
832  Int_t isPhiKKpi=isSel&4;
833  Int_t isPhipiKK=isSel&8;
834  Int_t isK0starKKpi=isSel&16;
835  Int_t isK0starpiKK=isSel&32;
836 
837 
838  if(fDsChannel==kPhi && (isPhiKKpi==0 && isPhipiKK==0)) return;
839  if(fDsChannel==kK0star && (isK0starKKpi==0 && isK0starpiKK==0)) return;
840 
841  if (optDecay==1){
842  if(isKKpi && fPartOrAndAntiPart*d->GetCharge()>=0) {
843  if(fDsChannel==kPhi && isPhiKKpi==0) return;
844  if(fDsChannel==kK0star && isK0starKKpi==0) return;
845 
846  fMassHist[index]->Fill(masses[0]);
847 
848  if(fReadMC){
849  if(labDs>=0){
850  Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
851  AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
852  Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
853  if(pdgCode0==321) {
854  fSigHist[index]->Fill(masses[0]); //signal
855  }else{
856  fRflHist[index]->Fill(masses[0]); //Reflected signal
857  }
858  }else{
859  fBkgHist[index]->Fill(masses[0]); // Background
860  }
861  }
862  }
863  }
864 
865  if (optDecay==0){
866  if(ispiKK && fPartOrAndAntiPart*d->GetCharge()>=0){
867  if(fDsChannel==kPhi && isPhipiKK==0) return;
868  if(fDsChannel==kK0star && isK0starpiKK==0) return;
869 
870  fMassHist[index]->Fill(masses[1]);
871 
872  if(fReadMC){
873  if(labDs>=0){
874  Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
875  AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
876  Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
877  if(pdgCode0==211) {
878  fSigHist[index]->Fill(masses[1]);
879  }else{
880  fRflHist[index]->Fill(masses[1]);
881  }
882  }else{
883  fBkgHist[index]->Fill(masses[1]);
884  }
885  }
886  }
887  }
888 }
889 
890 void AliAnalysisTaskSESignificance::FillD04p(AliAODRecoDecayHF* /*d*/,TClonesArray */*arrayMC*/,Int_t /*index*/,Int_t /*matchtoMC*/){
891  //D0->Kpipipi channel
892  AliInfo("D0 in 4 prongs channel not implemented\n");
893 
894 }
895 
896 void AliAnalysisTaskSESignificance::FillLambdac(AliAODRecoDecayHF* d,TClonesArray* arrayMC,Int_t index,Int_t isSel){
897 
898  // Mass hypothesis
899  Double_t masses[2];
900  Int_t pdgdaughtersLc[3]={2212,321,211}; //p,K,pi
901  masses[0]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersLc);
902  Int_t pdgdaughtersLc2[3]={211,321,2212}; //pi,K,p
903  masses[1]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersLc2);
904 
905 
906  if(fPartOrAndAntiPart==0 || fPartOrAndAntiPart==d->GetCharge()) {
907 
908  // isSel=1 : p K pi ; isSel=2 : pi K p ;
909  if(isSel==1 || isSel==3) fMassHist[index]->Fill(masses[0]);
910  if(isSel>=2) fMassHist[index]->Fill(masses[1]);
911 
912  // Check the MC truth
913  if(fReadMC){
914 
915  Int_t ispKpi = 0;
916  Int_t ispiKp = 0;
917  ispKpi = isSel&1;
918  ispiKp = isSel&2;
919  Int_t matchtoMC = -1;
920  //
921  Int_t pPDG = 2212; // p
922  Int_t kPDG = 321; // K
923  Int_t piPDG = 211; // pi
924  Int_t absPdgMom = 4122;
925  matchtoMC = d->MatchToMC(absPdgMom,arrayMC,fNProngs,pdgdaughtersLc);
926 
927  if(matchtoMC>=0){
928 
929  AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
930  Int_t pdgMC = dMC->GetPdgCode();
931  if (TMath::Abs(pdgMC)!=absPdgMom) AliInfo("What's up, isn't it a lambdac ?!");
932  Int_t labDau0 = ((AliAODTrack*)d->GetDaughter(0))->GetLabel();
933  AliAODMCParticle* p0 = (AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
934  Int_t pdgCode0 = TMath::Abs(p0->GetPdgCode());
935  Int_t labDau1 = ((AliAODTrack*)d->GetDaughter(1))->GetLabel();
936  AliAODMCParticle* p1 = (AliAODMCParticle*)arrayMC->UncheckedAt(labDau1);
937  Int_t pdgCode1 = TMath::Abs(p1->GetPdgCode());
938  Int_t labDau2 = ((AliAODTrack*)d->GetDaughter(2))->GetLabel();
939  AliAODMCParticle* p2 = (AliAODMCParticle*)arrayMC->UncheckedAt(labDau2);
940  Int_t pdgCode2 = TMath::Abs(p2->GetPdgCode());
941 
942  // Fill in the histograms in case of p K pi decays
943  if(ispKpi==1){
944  if(pdgCode0==pPDG && pdgCode1==kPDG && pdgCode2==piPDG){
945  fSigHist[index]->Fill(masses[0]);
946  } else {
947  fRflHist[index]->Fill(masses[0]);
948  }
949  }
950  // Fill in the histograms in case of pi K p decays
951  if(ispiKp==2){
952  if(pdgCode0==piPDG && pdgCode1==kPDG && pdgCode2==pPDG){
953  fSigHist[index]->Fill(masses[1]);
954  } else {
955  fRflHist[index]->Fill(masses[1]);
956  }
957  }
958  } else {
959  if(ispKpi==1) fBkgHist[index]->Fill(masses[0]);
960  if(ispiKp==2) fBkgHist[index]->Fill(masses[1]);
961  }
962  }
963  }
964 
965 
966 }
967 
968 
969 //________________________________________________________________________
971 {
972  // Terminate analysis
973  //
974  if(fDebug > 1) printf("AnalysisTaskSESignificance: Terminate() \n");
975 
976 
977  fOutput = dynamic_cast<TList*> (GetOutputData(1));
978  if (!fOutput) {
979  printf("ERROR: fOutput not available\n");
980  return;
981  }
982 
983  fCutList = dynamic_cast<TList*> (GetOutputData(2));
984  if (!fCutList) {
985  printf("ERROR: fCutList not available\n");
986  return;
987  }
988 
989  AliMultiDimVector* mdvtmp=(AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0");
990  if (!mdvtmp){
991  cout<<"multidimvec not found in TList"<<endl;
992  fCutList->ls();
993  return;
994  }
995  Int_t nHist=mdvtmp->GetNTotCells();
996  TCanvas *c1=new TCanvas("c1","Invariant mass distribution - loose cuts",500,500);
997  Bool_t drawn=kFALSE;
998  for(Int_t i=0;i<nHist;i++){
999 
1000  TString hisname;
1001  TString signame;
1002  TString bkgname;
1003  TString rflname;
1004 
1005  hisname.Form("hMass_%d",i);
1006  signame.Form("hSig_%d",i);
1007  bkgname.Form("hBkg_%d",i);
1008  rflname.Form("hRfl_%d",i);
1009 
1010  fMassHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
1011 
1012  if (!drawn && fMassHist[i]->GetEntries() > 0){
1013  c1->cd();
1014  fMassHist[i]->Draw();
1015  drawn=kTRUE;
1016  }
1017 
1018  if(fReadMC){
1019  fSigHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(signame.Data()));
1020  fBkgHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(bkgname.Data()));
1021  if(fDecChannel != AliAnalysisTaskSESignificance::kDplustoKpipi) fRflHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(rflname.Data()));
1022  }
1023 
1024  }
1025 
1026 
1027 
1028 
1029  return;
1030 }
1031 //_________________________________________________________________________________________________
1032 Int_t AliAnalysisTaskSESignificance::CheckOrigin(const AliAODMCParticle* mcPart, const TClonesArray* mcArray)const{
1033 
1034  //
1035  // checking whether the very mother of the D0 is a charm or a bottom quark
1036  //
1037 
1038  Int_t pdgGranma = 0;
1039  Int_t mother = 0;
1040  mother = mcPart->GetMother();
1041  Int_t istep = 0;
1042  while (mother >0 ){
1043  istep++;
1044  AliDebug(2,Form("mother at step %d = %d", istep, mother));
1045  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
1046  if(!mcGranma) break;
1047  pdgGranma = mcGranma->GetPdgCode();
1048  AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1049  Int_t abspdgGranma = TMath::Abs(pdgGranma);
1050  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1051  break;
1052  }
1053  mother = mcGranma->GetMother();
1054  }
1055  return pdgGranma;
1056 }
1057 //_________________________________________________________________________________________________
Float_t fLowmasslimit
upper inv mass limit for histos
virtual void GetCutVarsForOpt(AliAODRecoDecayHF *d, Float_t *vars, Int_t nvars, Int_t *pdgdaughters)=0
Int_t pdg
Int_t GetNVars() const
Definition: AliRDHFCuts.h:249
double Double_t
Definition: External.C:58
Float_t fVars[kMaxCutVar]
number of selection variables
AliRDHFCuts * fRDCuts
lower inv mass limit for histos
virtual void UserExec(Option_t *option)
const char * title
Definition: MakeQAPdf.C:27
Double_t DeltaInvMass() const
Bool_t HasSelectionBit(Int_t i) const
Int_t fNBins
array with values of cut variables
Int_t MatchToMC(Int_t pdgabs, Int_t pdgabs2prong, Int_t *pdgDg, Int_t *pdgDg2prong, TClonesArray *mcArray, Bool_t isV0=kFALSE) const
Int_t fPDGD0ToKpi[2]
PDG codes for the particles in the D* -> pi + D0 decay.
Double_t mass
Int_t fPDGdaughters[4]
number of prong of the decay channel
static Int_t CheckMatchingAODdeltaAODevents()
TH1F * fBkgHist[kMaxNHist]
!hist. for inv mass (bkg from MC truth)
FeedDownEnum fBFeedDown
flag to use selection bit (speed up candidates selection)
Int_t GetWhyRejection() const
Definition: AliRDHFCuts.h:313
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
TList * fOutput
! list send on output slot 0
TString * GetVarNames() const
Definition: AliRDHFCuts.h:250
Int_t fPartOrAndAntiPart
number of bins in the mass histograms
void FillDs(AliAODRecoDecayHF *d, TClonesArray *arrayMC, Int_t index, Int_t isSel, Int_t optDecay)
Int_t fSelectionlevel
AOD branch name for channel.
Bool_t FillRecoCasc(AliVEvent *event, AliAODRecoCascadeHF *rc, Bool_t isDStar, Bool_t recoSecVtx=kFALSE)
Bool_t fUseSelBit
flag for access to MC
Class for cuts on AOD reconstructed D+->Kpipi.
void FillLambdac(AliAODRecoDecayHF *d, TClonesArray *arrayMC, Int_t index, Int_t isSel)
int Int_t
Definition: External.C:63
Bool_t * GetIsUpperCut() const
Definition: AliRDHFCuts.h:260
TH1F * fRflHist[kMaxNHist]
!hist. for inv mass (bkg from MC truth)
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
Bool_t * GetVarsForOpt() const
Definition: AliRDHFCuts.h:251
Int_t CheckOrigin(const AliAODMCParticle *mcPart, const TClonesArray *mcArray) const
virtual void UserCreateOutputObjects()
Implementation of interface methods.
Int_t GetNVarsForOpt() const
Definition: AliRDHFCuts.h:252
ULong64_t * GetGlobalAddressesAboveCuts(const Float_t *values, Float_t pt, Int_t &nVals) const
Int_t fPDGmother
decay channel identifier
void FillD04p(AliAODRecoDecayHF *d, TClonesArray *arrayMC, Int_t index, Int_t isSel)
Int_t fPDGDStarToD0pi[2]
Ds resonant channel selected.
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
TString fBranchName
PDG codes of daughters.
void FillDplus(AliAODRecoDecayHF *d, TClonesArray *arrayMC, Int_t index, Int_t isSel)
Bool_t IsEventSelected(AliVEvent *event)
ULong64_t GetNTotCells() const
virtual void PrintAll() const
void FillD02p(AliAODRecoDecayHF *d, TClonesArray *arrayMC, Int_t index, Int_t isSel)
void SetMassLimits(Float_t range, Int_t pdg)
Bool_t IsSelected(TObject *obj)
Definition: AliRDHFCuts.h:290
virtual void Terminate(Option_t *option)
Int_t fNVars
selection level: kALL,kTracks,kCandidate
const char Option_t
Definition: External.C:48
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:248
bool Bool_t
Definition: External.C:53
AliAODRecoDecayHF2Prong * Get2Prong() const
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:310
Int_t PtBin(Double_t pt) const
TH1F * fMassHist[kMaxNHist]
Multidimvector container.
Int_t fDecChannel
flag to search for D from B decays
void FillDstar(AliAODRecoCascadeHF *dstarD0pi, TClonesArray *arrayMC, Int_t index, Int_t isSel)
Int_t fDsChannel
fill histograms with particle only (+1), antiparticle only (-1), both (0)
TH1F * fSigHist[kMaxNHist]
!hist. for inv mass (sig from MC truth)