AliPhysics  d497547 (d497547)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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),
130  fPartOrAndAntiPart(0),
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 
538  for (Int_t iProng = 0; iProng < nProng; iProng++) {
539  fHistNEvents->Fill(6);
540  d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iProng);
541 
542  AliAODRecoCascadeHF* DStarToD0pi = NULL;
543  AliAODRecoDecayHF2Prong* D0Particle = NULL;
544  fPDGDStarToD0pi[0] = 421; fPDGDStarToD0pi[1] = 211;
545  fPDGD0ToKpi[0] = 321; fPDGD0ToKpi[1] = 211;
546 
547  Bool_t isSelBit=kTRUE;
548  if(fUseSelBit){
549  if(fDecChannel==0) {
551  }else{
552  if(fDecChannel==1) {
554  }else{
555  if(fDecChannel==2) {
557  }else{
558  if(fDecChannel==3) {
560  }else{
562  }
563  }
564  }
565  }
566  }
567  if(!isSelBit) continue;
568 
569  if (fDecChannel==2) {
570  DStarToD0pi = (AliAODRecoCascadeHF*)arrayProng->At(iProng);
571  if (!DStarToD0pi->GetSecondaryVtx()) continue;
572  D0Particle = (AliAODRecoDecayHF2Prong*)DStarToD0pi->Get2Prong();
573  if (!D0Particle) continue;
574  }
575 
576  Bool_t isFidAcc = fRDCuts->IsInFiducialAcceptance(d->Pt(),d->Y(fPDGmother));
577  Int_t isSelected=fRDCuts->IsSelected(d,fSelectionlevel,aod);
578 
579  if(fReadMC && fBFeedDown!=kBoth && isSelected){
580  Int_t labD = d->MatchToMC(fPDGmother,arrayMC,fNProngs,fPDGdaughters);
581  if(labD>=0){
582  AliAODMCParticle *partD = (AliAODMCParticle*)arrayMC->At(labD);
583  Int_t pdgGranma = CheckOrigin(partD, arrayMC);
584  Int_t abspdgGranma = TMath::Abs(pdgGranma);
585  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
586  //feed down particle
587  AliDebug(2,Form("Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
588  fHistNEvents->Fill(7);
589  if(fBFeedDown==kCharmOnly) isSelected=kFALSE; //from beauty
590  }
591  else {
592  //prompt particle
593  fHistNEvents->Fill(6);
594  if(fBFeedDown==kBeautyOnly)isSelected=kFALSE;
595  }
596  }
597  }
598 
599  if(isSelected&&isFidAcc) {
600  fHistNEvents->Fill(2); // count selected with loosest cuts
601  if(fDebug>1) printf("+++++++Is Selected\n");
602 
603  Int_t nVals=0;
606  Int_t ptbin=fRDCuts->PtBin(d->Pt());
607  if(ptbin==-1) continue;
608  TString mdvname=Form("multiDimVectorPtBin%d",ptbin);
609  AliMultiDimVector* muvec=(AliMultiDimVector*)fCutList->FindObject(mdvname.Data());
610 
611  ULong64_t *addresses = muvec->GetGlobalAddressesAboveCuts(fVars,(Float_t)d->Pt(),nVals);
612  if(fDebug>1)printf("nvals = %d\n",nVals);
613  for(Int_t ivals=0;ivals<nVals;ivals++){
614  if(addresses[ivals]>=muvec->GetNTotCells()){
615  if (fDebug>1) printf("Overflow!!\n");
616  delete [] addresses;
617  return;
618  }
619 
620  fHistNEvents->Fill(3);
621 
622  //fill the histograms with the appropriate method
623  switch (fDecChannel){
624  case 0:
625  FillDplus(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
626  break;
627  case 1:
628  FillD02p(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
629  break;
630  case 2:
631  FillDstar(DStarToD0pi,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
632  break;
633  case 3:
634  if(isSelected&1){
635  FillDs(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected,1);
636  }
637  break;
638  case 4:
639  FillD04p(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
640  break;
641  case 5:
642  FillLambdac(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected);
643  break;
644  default:
645  break;
646  }
647 
648  }
649 
650  if (fDecChannel==3 && isSelected&2){
652  nVals=0;
654  delete [] addresses;
655  addresses = muvec->GetGlobalAddressesAboveCuts(fVars,(Float_t)d->Pt(),nVals);
656  if(fDebug>1)printf("nvals = %d\n",nVals);
657  for(Int_t ivals=0;ivals<nVals;ivals++){
658  if(addresses[ivals]>=muvec->GetNTotCells()){
659  if (fDebug>1) printf("Overflow!!\n");
660  delete [] addresses;
661  return;
662  }
663  FillDs(d,arrayMC,(Int_t)(ptbin*nHistpermv+addresses[ivals]),isSelected,0);
664 
665  }
666 
667  }
668 
669  delete [] addresses;
670  }// end if selected
671 
672  }
673 
674 
675  PostData(1,fOutput);
676  return;
677 }
678 
679 //***************************************************************************
680 
681 // Methods used in the UserExec
682 
683 
684 //********************************************************************************************
685 
686 //Methods to fill the istograms with MC information, one for each candidate
687 //NB: the implementation for each candidate is responsibility of the corresponding developer
688 
689 void AliAnalysisTaskSESignificance::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel){
690  //D+ channel
691  if(!isSel){
692  AliError("Candidate not selected\n");
693  return;
694  }
695 
696  if(fPartOrAndAntiPart*d->GetCharge()<0)return;
697  Int_t pdgdaughters[3] = {211,321,211};
698  Double_t mass=d->InvMass(3,(UInt_t*)pdgdaughters);
699 
700  fMassHist[index]->Fill(mass);
701 
702 
703  if(fReadMC){
704  Int_t lab=-1;
705  lab = d->MatchToMC(411,arrayMC,3,pdgdaughters);
706  if(lab>=0){ //signal
707  fSigHist[index]->Fill(mass);
708  } else{ //background
709  fBkgHist[index]->Fill(mass);
710  }
711  }
712 }
713 
714 
715 void AliAnalysisTaskSESignificance::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel){
716 
717  //D0->Kpi channel
718 
719  //mass histograms
720  Int_t pdgdaughtersD0[2]={211,321};//pi,K
721  Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
722  Double_t masses[2];
723  masses[0]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
724  masses[1]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
725 
726  if((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0) fMassHist[index]->Fill(masses[0]);
727  if(isSel>=2 && fPartOrAndAntiPart<=0) fMassHist[index]->Fill(masses[1]);
728 
729 
730 
731  //MC histograms
732  if(fReadMC){
733 
734  Int_t matchtoMC=-1;
735 
736  //D0
737  Int_t pdgdaughters[2];
738  pdgdaughters[0]=211;//pi
739  pdgdaughters[1]=321;//K
740 
741  matchtoMC = d->MatchToMC(fPDGmother,arrayMC,fNProngs,pdgdaughters);
742 
743  Int_t prongPdgPlus=421,prongPdgMinus=(-1)*421;
744  if((isSel==1 || isSel==3) && fPartOrAndAntiPart>=0){ //D0
745  if(matchtoMC>=0){
746  AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
747  Int_t pdgMC = dMC->GetPdgCode();
748 
749  if(pdgMC==prongPdgPlus) fSigHist[index]->Fill(masses[0]);
750  else fRflHist[index]->Fill(masses[0]);
751 
752  } else fBkgHist[index]->Fill(masses[0]);
753 
754  }
755  if(isSel>=2 && fPartOrAndAntiPart<=0){ //D0bar
756  if(matchtoMC>=0){
757  AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
758  Int_t pdgMC = dMC->GetPdgCode();
759 
760  if(pdgMC==prongPdgMinus) fSigHist[index]->Fill(masses[1]);
761  else fRflHist[index]->Fill(masses[1]);
762  } else fBkgHist[index]->Fill(masses[1]);
763  }
764  }
765 }
766 
767 void AliAnalysisTaskSESignificance::FillDstar(AliAODRecoCascadeHF* dstarD0pi,TClonesArray *arrayMC,Int_t index,Int_t isSel){
768  //D* channel
769 
770  AliInfo("Dstar selected\n");
771 
772  Double_t mass = dstarD0pi->DeltaInvMass();
773 
774  if((isSel>0) && TMath::Abs(fPartOrAndAntiPart)>=0) fMassHist[index]->Fill(mass);
775 
776  if(fReadMC) {
777  Int_t matchtoMC = -1;
778  matchtoMC = dstarD0pi->MatchToMC(413,421,fPDGDStarToD0pi, fPDGD0ToKpi, arrayMC);
779 
780  Int_t prongPdgDStarPlus=413;//,prongPdgDStarMinus=(-1)*413;
781 
782  if ((isSel>1) && TMath::Abs(fPartOrAndAntiPart)>=0) {
783  //D*+
784  if(matchtoMC>=0) {
785  AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
786  Int_t pdgMC = dMC->GetPdgCode();
787 
788  if (pdgMC==prongPdgDStarPlus) fSigHist[index]->Fill(mass);
789  else {
790  dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
791  mass = dstarD0pi->DeltaInvMass();
792  fRflHist[index]->Fill(mass);
793  dstarD0pi->SetCharge(-1*dstarD0pi->GetCharge());
794  }
795  }
796  else fBkgHist[index]->Fill(mass);
797  }
798  }
799 }
800 
801 
802 void AliAnalysisTaskSESignificance::FillDs(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t index,Int_t isSel,Int_t optDecay){
803 
804  // Fill Ds histos
805 
806 
807  Int_t pdgDsKKpi[3]={321,321,211};//K,K,pi
808  Int_t pdgDspiKK[3]={211,321,321};//pi,K,K
809  Double_t masses[2];
810  masses[0]=d->InvMass(fNProngs,(UInt_t*)pdgDsKKpi); //Ds
811  masses[1]=d->InvMass(fNProngs,(UInt_t*)pdgDspiKK); //Dsbar
812 
813  Int_t labDs=-1;
814  if(fReadMC){
815  labDs = d->MatchToMC(431,arrayMC,3,pdgDsKKpi);
816  }
817 
818  Int_t isKKpi=isSel&1;
819  Int_t ispiKK=isSel&2;
820  Int_t isPhiKKpi=isSel&4;
821  Int_t isPhipiKK=isSel&8;
822  Int_t isK0starKKpi=isSel&16;
823  Int_t isK0starpiKK=isSel&32;
824 
825 
826  if(fDsChannel==kPhi && (isPhiKKpi==0 && isPhipiKK==0)) return;
827  if(fDsChannel==kK0star && (isK0starKKpi==0 && isK0starpiKK==0)) return;
828 
829  if (optDecay==1){
830  if(isKKpi && fPartOrAndAntiPart*d->GetCharge()>=0) {
831  if(fDsChannel==kPhi && isPhiKKpi==0) return;
832  if(fDsChannel==kK0star && isK0starKKpi==0) return;
833 
834  fMassHist[index]->Fill(masses[0]);
835 
836  if(fReadMC){
837  if(labDs>=0){
838  Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
839  AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
840  Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
841  if(pdgCode0==321) {
842  fSigHist[index]->Fill(masses[0]); //signal
843  }else{
844  fRflHist[index]->Fill(masses[0]); //Reflected signal
845  }
846  }else{
847  fBkgHist[index]->Fill(masses[0]); // Background
848  }
849  }
850  }
851  }
852 
853  if (optDecay==0){
854  if(ispiKK && fPartOrAndAntiPart*d->GetCharge()>=0){
855  if(fDsChannel==kPhi && isPhipiKK==0) return;
856  if(fDsChannel==kK0star && isK0starpiKK==0) return;
857 
858  fMassHist[index]->Fill(masses[1]);
859 
860  if(fReadMC){
861  if(labDs>=0){
862  Int_t labDau0=((AliAODTrack*)d->GetDaughter(0))->GetLabel();
863  AliAODMCParticle* p=(AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
864  Int_t pdgCode0=TMath::Abs(p->GetPdgCode());
865  if(pdgCode0==211) {
866  fSigHist[index]->Fill(masses[1]);
867  }else{
868  fRflHist[index]->Fill(masses[1]);
869  }
870  }else{
871  fBkgHist[index]->Fill(masses[1]);
872  }
873  }
874  }
875  }
876 }
877 
878 void AliAnalysisTaskSESignificance::FillD04p(AliAODRecoDecayHF* /*d*/,TClonesArray */*arrayMC*/,Int_t /*index*/,Int_t /*matchtoMC*/){
879  //D0->Kpipipi channel
880  AliInfo("D0 in 4 prongs channel not implemented\n");
881 
882 }
883 
884 void AliAnalysisTaskSESignificance::FillLambdac(AliAODRecoDecayHF* d,TClonesArray* arrayMC,Int_t index,Int_t isSel){
885 
886  // Mass hypothesis
887  Double_t masses[2];
888  Int_t pdgdaughtersLc[3]={2212,321,211}; //p,K,pi
889  masses[0]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersLc);
890  Int_t pdgdaughtersLc2[3]={211,321,2212}; //pi,K,p
891  masses[1]=d->InvMass(fNProngs,(UInt_t*)pdgdaughtersLc2);
892 
893 
894  if(fPartOrAndAntiPart==0 || fPartOrAndAntiPart==d->GetCharge()) {
895 
896  // isSel=1 : p K pi ; isSel=2 : pi K p ;
897  if(isSel==1 || isSel==3) fMassHist[index]->Fill(masses[0]);
898  if(isSel>=2) fMassHist[index]->Fill(masses[1]);
899 
900  // Check the MC truth
901  if(fReadMC){
902 
903  Int_t ispKpi = 0;
904  Int_t ispiKp = 0;
905  ispKpi = isSel&1;
906  ispiKp = isSel&2;
907  Int_t matchtoMC = -1;
908  //
909  Int_t pPDG = 2212; // p
910  Int_t kPDG = 321; // K
911  Int_t piPDG = 211; // pi
912  Int_t absPdgMom = 4122;
913  matchtoMC = d->MatchToMC(absPdgMom,arrayMC,fNProngs,pdgdaughtersLc);
914 
915  if(matchtoMC>=0){
916 
917  AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
918  Int_t pdgMC = dMC->GetPdgCode();
919  if (TMath::Abs(pdgMC)!=absPdgMom) AliInfo("What's up, isn't it a lambdac ?!");
920  Int_t labDau0 = ((AliAODTrack*)d->GetDaughter(0))->GetLabel();
921  AliAODMCParticle* p0 = (AliAODMCParticle*)arrayMC->UncheckedAt(labDau0);
922  Int_t pdgCode0 = TMath::Abs(p0->GetPdgCode());
923  Int_t labDau1 = ((AliAODTrack*)d->GetDaughter(1))->GetLabel();
924  AliAODMCParticle* p1 = (AliAODMCParticle*)arrayMC->UncheckedAt(labDau1);
925  Int_t pdgCode1 = TMath::Abs(p1->GetPdgCode());
926  Int_t labDau2 = ((AliAODTrack*)d->GetDaughter(2))->GetLabel();
927  AliAODMCParticle* p2 = (AliAODMCParticle*)arrayMC->UncheckedAt(labDau2);
928  Int_t pdgCode2 = TMath::Abs(p2->GetPdgCode());
929 
930  // Fill in the histograms in case of p K pi decays
931  if(ispKpi==1){
932  if(pdgCode0==pPDG && pdgCode1==kPDG && pdgCode2==piPDG){
933  fSigHist[index]->Fill(masses[0]);
934  } else {
935  fRflHist[index]->Fill(masses[0]);
936  }
937  }
938  // Fill in the histograms in case of pi K p decays
939  if(ispiKp==2){
940  if(pdgCode0==piPDG && pdgCode1==kPDG && pdgCode2==pPDG){
941  fSigHist[index]->Fill(masses[1]);
942  } else {
943  fRflHist[index]->Fill(masses[1]);
944  }
945  }
946  } else {
947  if(ispKpi==1) fBkgHist[index]->Fill(masses[0]);
948  if(ispiKp==2) fBkgHist[index]->Fill(masses[1]);
949  }
950  }
951  }
952 
953 
954 }
955 
956 
957 //________________________________________________________________________
959 {
960  // Terminate analysis
961  //
962  if(fDebug > 1) printf("AnalysisTaskSESignificance: Terminate() \n");
963 
964 
965  fOutput = dynamic_cast<TList*> (GetOutputData(1));
966  if (!fOutput) {
967  printf("ERROR: fOutput not available\n");
968  return;
969  }
970 
971  fCutList = dynamic_cast<TList*> (GetOutputData(2));
972  if (!fCutList) {
973  printf("ERROR: fCutList not available\n");
974  return;
975  }
976 
977  AliMultiDimVector* mdvtmp=(AliMultiDimVector*)fCutList->FindObject("multiDimVectorPtBin0");
978  if (!mdvtmp){
979  cout<<"multidimvec not found in TList"<<endl;
980  fCutList->ls();
981  return;
982  }
983  Int_t nHist=mdvtmp->GetNTotCells();
984  TCanvas *c1=new TCanvas("c1","Invariant mass distribution - loose cuts",500,500);
985  Bool_t drawn=kFALSE;
986  for(Int_t i=0;i<nHist;i++){
987 
988  TString hisname;
989  TString signame;
990  TString bkgname;
991  TString rflname;
992 
993  hisname.Form("hMass_%d",i);
994  signame.Form("hSig_%d",i);
995  bkgname.Form("hBkg_%d",i);
996  rflname.Form("hRfl_%d",i);
997 
998  fMassHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(hisname.Data()));
999 
1000  if (!drawn && fMassHist[i]->GetEntries() > 0){
1001  c1->cd();
1002  fMassHist[i]->Draw();
1003  drawn=kTRUE;
1004  }
1005 
1006  if(fReadMC){
1007  fSigHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(signame.Data()));
1008  fBkgHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(bkgname.Data()));
1009  if(fDecChannel != AliAnalysisTaskSESignificance::kDplustoKpipi) fRflHist[i]=dynamic_cast<TH1F*>(fOutput->FindObject(rflname.Data()));
1010  }
1011 
1012  }
1013 
1014 
1015 
1016 
1017  return;
1018 }
1019 //_________________________________________________________________________________________________
1020 Int_t AliAnalysisTaskSESignificance::CheckOrigin(const AliAODMCParticle* mcPart, const TClonesArray* mcArray)const{
1021 
1022  //
1023  // checking whether the very mother of the D0 is a charm or a bottom quark
1024  //
1025 
1026  Int_t pdgGranma = 0;
1027  Int_t mother = 0;
1028  mother = mcPart->GetMother();
1029  Int_t istep = 0;
1030  while (mother >0 ){
1031  istep++;
1032  AliDebug(2,Form("mother at step %d = %d", istep, mother));
1033  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
1034  if(!mcGranma) break;
1035  pdgGranma = mcGranma->GetPdgCode();
1036  AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1037  Int_t abspdgGranma = TMath::Abs(pdgGranma);
1038  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1039  break;
1040  }
1041  mother = mcGranma->GetMother();
1042  }
1043  return pdgGranma;
1044 }
1045 //_________________________________________________________________________________________________
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:240
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:304
TList * fOutput
! list send on output slot 0
TString * GetVarNames() const
Definition: AliRDHFCuts.h:241
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 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:251
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:242
Int_t CheckOrigin(const AliAODMCParticle *mcPart, const TClonesArray *mcArray) const
virtual void UserCreateOutputObjects()
Implementation of interface methods.
Int_t GetNVarsForOpt() const
Definition: AliRDHFCuts.h:243
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:281
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:239
bool Bool_t
Definition: External.C:53
AliAODRecoDecayHF2Prong * Get2Prong() const
virtual Bool_t IsInFiducialAcceptance(Double_t, Double_t) const
Definition: AliRDHFCuts.h:301
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)