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