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