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