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