]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/FlavourJetTasks/AliAnalysisTaskSEDmesonsFilterCJ.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
dont create any files in place to avoid problems with eos
[u/mrichter/AliRoot.git] / PWGJE / FlavourJetTasks / AliAnalysisTaskSEDmesonsFilterCJ.cxx
CommitLineData
c683985a 1/**************************************************************************
2* Copyright(c) 1998-2009, 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// Task for selecting D mesons to be used as an input for D-Jet correlations
18//
19//-----------------------------------------------------------------------
20// Authors:
21// C. Bianchin (Utrecht University) chiara.bianchin@cern.ch
22// A.Grelli (Utrecht University) a.grelli@uu.nl
23// Xiaoming Zhang (LBNL) XMZhang@lbl.gov
24//-----------------------------------------------------------------------
25
26#include <TDatabasePDG.h>
27#include <TParticle.h>
28#include <TVector3.h>
29#include "TROOT.h"
30#include <TH3F.h>
31
32#include "AliRDHFCutsDStartoKpipi.h"
33#include "AliRDHFCutsD0toKpi.h"
34#include "AliAODMCHeader.h"
35#include "AliAODHandler.h"
36#include "AliAnalysisManager.h"
37#include "AliLog.h"
38#include "AliAODVertex.h"
39#include "AliAODRecoDecay.h"
40#include "AliAODRecoCascadeHF.h"
41#include "AliAODRecoDecayHF2Prong.h"
42#include "AliESDtrack.h"
43#include "AliAODMCParticle.h"
44#include "AliAnalysisTaskSEDmesonsFilterCJ.h"
45
46ClassImp(AliAnalysisTaskSEDmesonsFilterCJ)
47
48//_______________________________________________________________________________
49
50AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ() :
51AliAnalysisTaskSE(),
52fUseMCInfo(kFALSE),
53fUseReco(kTRUE),
54fCandidateType(0),
55fCandidateName(""),
56fPDGmother(0),
57fNProngs(0),
58fBranchName(""),
59fOutput(0),
60fCuts(0),
61fMinMass(0.),
62fMaxMass(0.),
bbb94467 63fCandidateArray(0),
64fSideBandArray(0)
c683985a 65
66{
67 //
68 // Default constructor
69 //
70
71 for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
72 for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
73}
74
75//_______________________________________________________________________________
76
77AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ(const char *name, AliRDHFCuts *cuts, ECandidateType candtype) :
78AliAnalysisTaskSE(name),
79fUseMCInfo(kFALSE),
80fUseReco(kTRUE),
81fCandidateType(candtype),
82fCandidateName(""),
83fPDGmother(0),
84fNProngs(0),
85fBranchName(""),
86fOutput(0),
87fCuts(cuts),
88fMinMass(0.),
89fMaxMass(0.),
bbb94467 90fCandidateArray(0),
91fSideBandArray(0)
c683985a 92{
93 //
94 // Constructor. Initialization of Inputs and Outputs
95 //
96
97 Info("AliAnalysisTaskSEDmesonsFilterCJ","Calling Constructor");
98
99 for (Int_t i=4; i--;) fPDGdaughters[i] = 0;
100 for (Int_t i=30; i--;) fSigmaD0[i] = 0.;
101
102 const Int_t nptbins = fCuts->GetNPtBins();
103 Float_t defaultSigmaD013[13] = { 0.012, 0.012, 0.012, 0.015, 0.015, 0.018, 0.018, 0.020, 0.020, 0.030, 0.030, 0.037, 0.040 };
104
105 switch (fCandidateType) {
106 case 0 :
107 fCandidateName = "D0";
108 fPDGmother = 421;
109 fNProngs = 2;
110 fPDGdaughters[0] = 211; // pi
111 fPDGdaughters[1] = 321; // K
112 fPDGdaughters[2] = 0; // empty
113 fPDGdaughters[3] = 0; // empty
114 fBranchName = "D0toKpi";
115 break;
116 case 1 :
117 fCandidateName = "Dstar";
118 fPDGmother = 413;
119 fNProngs = 3;
120 fPDGdaughters[1] = 211; // pi soft
121 fPDGdaughters[0] = 421; // D0
122 fPDGdaughters[2] = 211; // pi fromD0
123 fPDGdaughters[3] = 321; // K from D0
124 fBranchName = "Dstar";
125
126 if (nptbins<=13) {
127 for (Int_t ipt=0; ipt<nptbins;ipt++) fSigmaD0[ipt] = defaultSigmaD013[ipt];
128 } else {
129 AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
130 }
131 break;
132 default :
133 printf("%d not accepted!!\n",fCandidateType);
134 break;
135 }
136
137 if (fCandidateType==kD0toKpi) SetMassLimits(0.15, fPDGmother);
138 if (fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
139
140 DefineOutput(1, TList::Class()); // histos
141 DefineOutput(2, AliRDHFCuts::Class()); // my cuts
bbb94467 142 DefineOutput(3, TClonesArray::Class()); //array of candidates
143 DefineOutput(4, TClonesArray::Class()); //array of SB candidates
c683985a 144}
145
146//_______________________________________________________________________________
147
148AliAnalysisTaskSEDmesonsFilterCJ::~AliAnalysisTaskSEDmesonsFilterCJ()
149{
150 //
151 // destructor
152 //
153
154 Info("~AliAnalysisTaskSEDmesonsFilterCJ","Calling Destructor");
155
156 if (fOutput) { delete fOutput; fOutput = 0; }
157 if (fCuts) { delete fCuts; fCuts = 0; }
158 if (fCandidateArray) { delete fCandidateArray; fCandidateArray = 0; }
bbb94467 159 delete fSideBandArray;
c683985a 160
161}
162
163//_______________________________________________________________________________
164
165void AliAnalysisTaskSEDmesonsFilterCJ::Init()
166{
167 //
168 // Initialization
169 //
170
171 if(fDebug>1) printf("AnalysisTaskSEDmesonsForJetCorrelations::Init() \n");
172
173 switch (fCandidateType) {
bbb94467 174 case 0:
175 {
c683985a 176 AliRDHFCutsD0toKpi* copyfCutsDzero = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
177 copyfCutsDzero->SetName("AnalysisCutsDzero");
178 PostData(2, copyfCutsDzero); // Post the data
bbb94467 179 } break;
180 case 1:
181 {
182 AliRDHFCutsDStartoKpipi* copyfCutsDstar = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
183 copyfCutsDstar->SetName("AnalysisCutsDStar");
184 PostData(2, copyfCutsDstar); // Post the cuts
185 } break;
186 default: return;
c683985a 187 }
188
189 return;
190}
191
192//_______________________________________________________________________________
193
194void AliAnalysisTaskSEDmesonsFilterCJ::UserCreateOutputObjects()
195{
196 //
197 // output
198 //
199
200 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
201
202 fOutput = new TList(); fOutput->SetOwner();
203 DefineHistoForAnalysis(); // define histograms
204
bbb94467 205 if (fCandidateType==kD0toKpi){
206 fCandidateArray = new TClonesArray("AliAODRecoDecayHF",0);
207 fSideBandArray = new TClonesArray("AliAODRecoDecayHF",0);
208 }
c683985a 209
bbb94467 210 if (fCandidateType==kDstartoKpipi) {
211 fCandidateArray = new TClonesArray("AliAODRecoCascadeHF",0);
212 fSideBandArray = new TClonesArray("AliAODRecoCascadeHF",0);
c683985a 213 }
214
bbb94467 215 fCandidateArray->SetOwner();
216 fCandidateArray->SetName(Form("fCandidateArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen"));
217
218 //this is used for the DStar side bands and MC!
219 fSideBandArray->SetOwner();
220 fSideBandArray->SetName(Form("fSideBandArray%s%s",fCandidateName.Data(),fUseReco ? "rec" : "gen"));
221
c683985a 222 PostData(1, fOutput);
bbb94467 223 PostData(3, fCandidateArray);
224 PostData(4, fSideBandArray);
225
c683985a 226 return;
227}
228
229//_______________________________________________________________________________
230
231void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *){
232 //
233 // user exec
234 //
c683985a 235 // Load the event
236 AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
237
238 TClonesArray *arrayDStartoD0pi = 0;
239 if (!aodEvent && AODEvent() && IsStandardAOD()) {
240
241 // In case there is an AOD handler writing a standard AOD, use the AOD
242 // event in memory rather than the input (ESD) event.
243 aodEvent = dynamic_cast<AliAODEvent*>(AODEvent());
244
245 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
246 // have to taken from the AOD event hold by the AliAODExtension
247 AliAODHandler *aodHandler = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
248 if(aodHandler->GetExtensions()) {
249 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
250 AliAODEvent *aodFromExt = ext->GetAOD();
251 arrayDStartoD0pi = (TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
252 }
253 } else {
254 arrayDStartoD0pi = (TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
255 }
256
257 if (!arrayDStartoD0pi) {
258 AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
259 return;
260 } else {
261 AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
262 }
263
264 TClonesArray* mcArray = 0x0;
265 if (fUseMCInfo) {
266 mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
267 if (!mcArray) {
268 printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
269 return;
270 }
271 }
272
273 //Histograms
274 TH1I* hstat = (TH1I*)fOutput->FindObject("hstat");
275 TH1F* hnSBCandEv=(TH1F*)fOutput->FindObject("hnSBCandEv");
276 TH1F* hnCandEv=(TH1F*)fOutput->FindObject("hnCandEv");
277 TH2F* hInvMassptD = (TH2F*)fOutput->FindObject("hInvMassptD");
278
279 TH1F* hPtPion=0x0;
280 if (fCandidateType==kDstartoKpipi) hPtPion = (TH1F*)fOutput->FindObject("hPtPion");
281 hstat->Fill(0);
282
283 // fix for temporary bug in ESDfilter
284 // the AODs with null vertex pointer didn't pass the PhysSel
285 if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
286
287 //Event selection
288 Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
289 //TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
290 if (!iseventselected) return;
291 hstat->Fill(1);
292
293 const Int_t nD = arrayDStartoD0pi->GetEntriesFast();
294 if(fUseReco) hstat->Fill(2, nD);
295
296 //D* and D0 prongs needed to MatchToMC method
297 Int_t pdgDgDStartoD0pi[2] = { 421, 211 }; // D0,pi
298 Int_t pdgDgD0toKpi[2] = { 321, 211 }; // K, pi
299
300 //D0 from D0 bar
301 Int_t pdgdaughtersD0[2] = { 211, 321 }; // pi,K
302 Int_t pdgdaughtersD0bar[2] = { 321, 211 }; // K,pi
303
304 Int_t iCand =0;
305 Int_t iSBCand=0;
306 Int_t isSelected = 0;
307 AliAODRecoDecayHF *charmCand = 0;
bbb94467 308 AliAODRecoCascadeHF *dstar = 0;
c683985a 309 AliAODMCParticle *charmPart = 0;
310 Bool_t isMCBkg=kFALSE;
311
312 Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
313
314 Int_t mcLabel = -9999;
315 Int_t pdgMeson = 413;
316 if (fCandidateType==kD0toKpi) pdgMeson = 421;
317
bbb94467 318 //clear the TClonesArray from the previous event
319 fCandidateArray->Clear();
320 fSideBandArray->Clear();
321
c683985a 322 for (Int_t icharm=0; icharm<nD; icharm++) { //loop over D candidates
323 charmCand = (AliAODRecoDecayHF*)arrayDStartoD0pi->At(icharm); // D candidates
324 if (!charmCand) continue;
325
8577f7bc 326 TString smcTruth="S";
327
bbb94467 328 if (fCandidateType==kDstartoKpipi) dstar = (AliAODRecoCascadeHF*)charmCand;
c683985a 329
330 if (fUseMCInfo) { // Look in MC, try to simulate the z
331 if (fCandidateType==kDstartoKpipi) {
bbb94467 332 mcLabel = dstar->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
c683985a 333 }
334
335 if (fCandidateType==kD0toKpi)
336 mcLabel = charmCand->MatchToMC(421,mcArray,fNProngs,fPDGdaughters);
337
338 if (mcLabel<=0) isMCBkg=kTRUE;
339 else hstat->Fill(2);
340 if (!isMCBkg) charmPart=(AliAODMCParticle*)mcArray->At(mcLabel);
8577f7bc 341
342 if (isMCBkg) smcTruth="B";
343
c683985a 344 }
345
346 Double_t ptD = charmCand->Pt();
347
348 // region of interest + cuts
349 if (!fCuts->IsInFiducialAcceptance(ptD,charmCand->Y(pdgMeson))) continue;
350
351 if(!fUseMCInfo && fCandidateType==kDstartoKpipi){
352 //select by track cuts the side band candidates (don't want mass cut)
353 isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kTracks,aodEvent);
354 if (!isSelected) continue;
355 //add a reasonable cut on the invariant mass (e.g. (+-2\sigma, +-10 \sigma), with \sigma = fSigmaD0[bin])
356 Int_t bin = fCuts->PtBin(ptD);
357 if(bin<0 || bin>=fCuts->GetNPtBins()) {
358 AliError(Form("Pt %.3f out of bounds",ptD));
359 continue;
360 }
c683985a 361 //if data and Dstar from D0 side band
bbb94467 362 if (((dstar->InvMassD0()<=(mPDGD0-3.*fSigmaD0[bin])) && (dstar->InvMassD0()>(mPDGD0-10.*fSigmaD0[bin]))) /*left side band*/|| ((dstar->InvMassD0()>=(mPDGD0+3.*fSigmaD0[bin])) && (dstar->InvMassD0()<(mPDGD0+10.*fSigmaD0[bin])))/*right side band*/){
c683985a 363
bbb94467 364 new ((*fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(*dstar);
c683985a 365 iSBCand++;
366 }
367 }
368 //candidate selected by cuts and PID
369 isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kAll,aodEvent); //selected
370 if (!isSelected) continue;
371
372 //for data and MC signal fill fCandidateArray
373 if(!fUseMCInfo || (fUseMCInfo && !isMCBkg)){
374 // for data or MC with the requirement fUseReco fill with candidates
375 if(fUseReco) {
bbb94467 376 if (fCandidateType==kDstartoKpipi){
377 new ((*fCandidateArray)[iCand]) AliAODRecoCascadeHF(*dstar);
378 AliInfo(Form("Dstar delta mass = %f",dstar->DeltaInvMass()));
379 } else{
380 new ((*fCandidateArray)[iCand]) AliAODRecoDecayHF(*charmCand);
381 //Printf("Filling reco");
382 }
383
c683985a 384 hstat->Fill(3);
385 }
386 // for MC with requirement particle level fill with AliAODMCParticle
387 else if (fUseMCInfo) {
388 new ((*fCandidateArray)[iCand]) AliAODMCParticle(*charmPart);
389 //Printf("Filling gen");
390 hstat->Fill(3);
391 }
392
393 iCand++;
394 }
395 //for MC background fill fSideBandArray (which is instead filled above for DStar in case of data for the side bands candidates)
396 else if(fUseReco){
bbb94467 397 if (fCandidateType==kDstartoKpipi){
398 new ((*fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(*dstar);
399 }
400 if (fCandidateType==kD0toKpi){
401 new ((*fSideBandArray)[iSBCand]) AliAODRecoDecayHF(*charmCand);
402 }
c683985a 403 iSBCand++;
404 }
405
406
407 Double_t masses[2];
408 if (fCandidateType==kDstartoKpipi) { //D*->D0pi->Kpipi
409
410 //softpion from D* decay
bbb94467 411 AliAODTrack *track2 = (AliAODTrack*)dstar->GetBachelor();
c683985a 412
413 // select D* in the D0 window.
414 // In the cut object window is loose to allow for side bands
415
416
417 // retrieve the corresponding pt bin for the candidate
418 // and set the expected D0 width (x3)
419 // static const Int_t n = fCuts->GetNPtBins();
420 Int_t bin = fCuts->PtBin(ptD);
421 if(bin<0 || bin>=fCuts->GetNPtBins()) {
422 AliError(Form("Pt %.3f out of bounds",ptD));
423 continue;
424 }
425
426 AliInfo(Form("Pt bin %d and sigma D0 %.4f",bin,fSigmaD0[bin]));
427 //consider the Dstar candidates only if the mass of the D0 is in 3 sigma wrt the PDG value
bbb94467 428 if ((dstar->InvMassD0()>=(mPDGD0-3.*fSigmaD0[bin])) && (dstar->InvMassD0()<=(mPDGD0+3.*fSigmaD0[bin]))) {
429 masses[0] = dstar->DeltaInvMass(); //D*
c683985a 430 masses[1] = 0.; //dummy for D*
431
432 //D* delta mass
433 hInvMassptD->Fill(masses[0], ptD); // 2 D slice for pt bins
434
435 // D* pt and soft pion pt for good candidates
436 Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
bbb94467 437 Double_t invmassDelta = dstar->DeltaInvMass();
c683985a 438 if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0021) hPtPion->Fill(track2->Pt());
439 }
8577f7bc 440
441 if (fUseMCInfo){ //fill histograms of kinematics, using MC truth
442 //get histos
443 TH2F *halphaDD = (TH2F*)fOutput->FindObject(Form("halphaDD%s",smcTruth.Data()));
444 TH2F *halphaDpis = (TH2F*)fOutput->FindObject(Form("halphaDpis%s",smcTruth.Data()));
445 TH2F *halphaDpi = (TH2F*)fOutput->FindObject(Form("halphaDpi%s",smcTruth.Data()));
446 TH2F *halphaDK = (TH2F*)fOutput->FindObject(Form("halphaDK%s",smcTruth.Data()));
447
448 TH2F *hdeltaRDD = (TH2F*)fOutput->FindObject(Form("hdeltaRDD%s",smcTruth.Data()));
449 TH2F *hdeltaRDpis = (TH2F*)fOutput->FindObject(Form("hdeltaRDpis%s",smcTruth.Data()));
450 TH2F *hdeltaRDpi = (TH2F*)fOutput->FindObject(Form("hdeltaRDpi%s",smcTruth.Data()));
451 TH2F *hdeltaRDK = (TH2F*)fOutput->FindObject(Form("hdeltaRDK%s",smcTruth.Data()));
452
453 Double_t aD = dstar->Phi(),
454 apis= track2->Phi();
455
456 AliAODRecoDecayHF2Prong* D0fromDstar=dstar->Get2Prong();
457 Double_t aD0 = D0fromDstar->Phi();
458 Int_t isD0= D0fromDstar->Charge()>0 ? kTRUE : kFALSE;
459 Double_t aK = isD0 ? D0fromDstar->PhiProng(0) : D0fromDstar->PhiProng(1),
460 api= isD0 ? D0fromDstar->PhiProng(1) : D0fromDstar->PhiProng(0);
461 Double_t dRDD0 = DeltaR(dstar,D0fromDstar),
462 dRDpis = DeltaR(dstar,track2),
463 dRDpi = DeltaR(dstar, isD0 ? (AliVParticle*)D0fromDstar->GetDaughter(1) : (AliVParticle*)D0fromDstar->GetDaughter(0)),
464 dRDK = DeltaR(dstar, isD0 ? (AliVParticle*)D0fromDstar->GetDaughter(0) : (AliVParticle*)D0fromDstar->GetDaughter(1));
465
466 halphaDD-> Fill(aD-aD0,ptD);
467 halphaDpis->Fill(aD-apis,ptD);
468 halphaDpi-> Fill(aD-api,ptD);
469 halphaDK-> Fill(aD-aK,ptD);
470
471 hdeltaRDD-> Fill(dRDD0,ptD);
472 hdeltaRDpis->Fill(dRDpis,ptD);
473 hdeltaRDpi-> Fill(dRDpi,ptD);
474 hdeltaRDK-> Fill(dRDK,ptD);
475
476 }
477
c683985a 478 } //Dstar specific
479
480 if (fCandidateType==kD0toKpi) { //D0->Kpi
481
482 //needed quantities
483 masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
484 masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
485 hstat->Fill(3);
486
487 // mass vs pt
488 if (isSelected==1 || isSelected==3) hInvMassptD->Fill(masses[0],ptD);
489 if (isSelected>=2) hInvMassptD->Fill(masses[1],ptD);
8577f7bc 490
491 if (fUseMCInfo) { //fill histograms of kinematics, using MC truth
492
493 Double_t aD = charmCand->Phi();
494 Double_t adaugh[2]={charmCand->PhiProng(0),charmCand->PhiProng(1)};
495 AliAODTrack* p0=(AliAODTrack*)charmCand->GetDaughter(0);
496 AliAODTrack* p1=(AliAODTrack*)charmCand->GetDaughter(1);
497 Float_t dR0 = DeltaR(charmCand, p0), dR1 = DeltaR(charmCand, p1);
498 Bool_t isD0=kFALSE;
499 if(mcLabel==421) isD0=kTRUE;
500 if(mcLabel==-421) isD0=kFALSE;
501
502 if(isMCBkg) { //background
503 TH2F *halphaDpi = (TH2F*)fOutput->FindObject(Form("halphaDpi%s",smcTruth.Data()));
504 TH2F *halphaDK = (TH2F*)fOutput->FindObject(Form("halphaDK%s",smcTruth.Data()));
505
506 TH2F *hdeltaRDpi = (TH2F*)fOutput->FindObject(Form("hdeltaRDpi%s",smcTruth.Data()));
507 TH2F *hdeltaRDK = (TH2F*)fOutput->FindObject(Form("hdeltaRDK%s",smcTruth.Data()));
508
509
510 if (isSelected==1 || isSelected==3) { // selected as D0
511 halphaDK->Fill(aD-adaugh[0],ptD);
512 halphaDpi->Fill(aD-adaugh[1],ptD);
513
514 hdeltaRDK->Fill(dR0,ptD);
515 hdeltaRDpi->Fill(dR1,ptD);
516
517 }
518 if (isSelected>=2) { //selected as D0bar
519 halphaDpi->Fill(aD-adaugh[0],ptD);
520 halphaDK->Fill(aD-adaugh[1],ptD);
521
522 hdeltaRDpi->Fill(dR0,ptD);
523 hdeltaRDK->Fill(dR1,ptD);
524
525 }
526
527 }else{ //signal and reflections
528 TH2F *halphaDpiS = (TH2F*)fOutput->FindObject("halphaDpiS");
529 TH2F *halphaDKS = (TH2F*)fOutput->FindObject("halphaDKS");
530 TH2F *halphaDpiR = (TH2F*)fOutput->FindObject("halphaDpiR");
531 TH2F *halphaDKR = (TH2F*)fOutput->FindObject("halphaDKR");
532
533 TH2F *hdeltaRDpiS = (TH2F*)fOutput->FindObject("hdeltaRDpiS");
534 TH2F *hdeltaRDKS = (TH2F*)fOutput->FindObject("hdeltaRDKS");
535 TH2F *hdeltaRDpiR = (TH2F*)fOutput->FindObject("hdeltaRDpiR");
536 TH2F *hdeltaRDKR = (TH2F*)fOutput->FindObject("hdeltaRDKR");
537
538 if(isD0) { //D0
539 halphaDKS->Fill(aD-adaugh[0],ptD);
540 halphaDpiS->Fill(aD-adaugh[1],ptD);
541
542 hdeltaRDKS->Fill(dR0,ptD);
543 hdeltaRDpiS->Fill(dR1,ptD);
544 if(isSelected>=2){ //selected as D0bar
545 halphaDpiR->Fill(aD-adaugh[0],ptD);
546 halphaDKR->Fill(aD-adaugh[1],ptD);
547
548 hdeltaRDpiR->Fill(dR0,ptD);
549 hdeltaRDKR->Fill(dR1,ptD);
550 }
551 } else { //D0bar
552 halphaDKS->Fill(aD-adaugh[1],ptD);
553 halphaDpiS->Fill(aD-adaugh[0],ptD);
554
555 hdeltaRDKS->Fill(dR1,ptD);
556 hdeltaRDpiS->Fill(dR0,ptD);
557
558 if(isSelected>=2){ //selected as D0bar
559 halphaDpiR->Fill(aD-adaugh[1],ptD);
560 halphaDKR->Fill(aD-adaugh[0],ptD);
561
562 hdeltaRDpiR->Fill(dR1,ptD);
563 hdeltaRDKR->Fill(dR0,ptD);
564 }
565 }
566
567 } //end signal and reflections
568
569
570 }// end MC
571
572
c683985a 573 } //D0 specific
574
575 charmCand = 0;
576 isMCBkg=kFALSE;
577 } // end of D cand loop
578
579 hnCandEv->Fill(fCandidateArray->GetEntriesFast());
580 if (fCandidateType==kDstartoKpipi || fUseMCInfo) {
581 Int_t nsbcand=fSideBandArray->GetEntriesFast();
582 hstat->Fill(4,nsbcand);
583 hnSBCandEv->Fill(nsbcand);
584 }
c3c1ab31 585 //Printf("N candidates selected %d, counter = %d",fCandidateArray->GetEntries(), iCand);
bbb94467 586 PostData(1, fOutput);
587 PostData(3, fCandidateArray);
588 PostData(4, fSideBandArray);
589
c683985a 590 return;
591}
592
593//_______________________________________________________________________________
594
595void AliAnalysisTaskSEDmesonsFilterCJ::Terminate(Option_t*)
596{
597 // The Terminate() function is the last function to be called during
598 // a query. It always runs on the client, it can be used to present
599 // the results graphically or save the results to file.
600
601 Info("Terminate"," terminate");
602 AliAnalysisTaskSE::Terminate();
603
604 fOutput = dynamic_cast<TList*>(GetOutputData(1));
605 if (!fOutput) {
606 printf("ERROR: fOutput not available\n");
607 return;
608 }
609
610 return;
611}
612
613//_______________________________________________________________________________
614
615void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t range, Int_t pdg)
616{
617 //
618 // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
619 //
620
621 Float_t mass = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdg))->Mass();
622
623 // compute the Delta mass for the D*
624 if (fCandidateType==kDstartoKpipi) mass -= TDatabasePDG::Instance()->GetParticle(421)->Mass();
625
626
627 fMinMass = mass - range;
628 fMaxMass = mass + range;
629
630 AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
631 if ((fMinMass<0.) || (fMaxMass<=0.) || (fMaxMass<fMinMass)) AliFatal("Wrong mass limits!\n");
632
633 return;
634}
635
636//_______________________________________________________________________________
637
638void AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits(Double_t lowlimit, Double_t uplimit)
639{
640 //
641 // AliAnalysisTaskSEDmesonsFilterCJ::SetMassLimits
642 //
643
644 if (uplimit>lowlimit) {
645 fMinMass = lowlimit;
646 fMaxMass = uplimit;
647 } else {
648 printf("Error! Lower limit larger than upper limit!\n");
649 fMinMass = uplimit - uplimit*0.2;
650 fMaxMass = uplimit;
651 }
652
653 return;
654}
655
656//_______________________________________________________________________________
657
658Bool_t AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar(Int_t nptbins, Float_t *width)
659{
660 //
661 // AliAnalysisTaskSEDmesonsFilterCJ::SetD0WidthForDStar
662 //
663
664 if (nptbins>30) {
665 AliInfo("Maximum number of bins allowed is 30!");
666 return kFALSE;
667 }
668
669 if (!width) return kFALSE;
670 for (Int_t ipt=0; ipt<nptbins; ipt++) fSigmaD0[ipt]=width[ipt];
671
672 return kTRUE;
673}
674
675//_______________________________________________________________________________
676
677Bool_t AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis()
678{
679 //
680 // AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis
681 //
682
683 // Statistics
684 TH1I* hstat = new TH1I("hstat","Statistics",5,-0.5,4.5);
685 hstat->GetXaxis()->SetBinLabel(1, "N ev anal");
686 hstat->GetXaxis()->SetBinLabel(2, "N ev sel");
687 if(fUseReco) hstat->GetXaxis()->SetBinLabel(3, "N cand");
688 else hstat->GetXaxis()->SetBinLabel(3, "N Gen D");
689 hstat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
690 if (fCandidateType==kDstartoKpipi) hstat->GetXaxis()->SetBinLabel(5, "N side band cand");
691 hstat->SetNdivisions(1);
692 fOutput->Add(hstat);
693
694 TH1F* hnCandEv=new TH1F("hnCandEv", "Number of candidates per event (after cuts);# cand/ev", 100, 0.,100.);
695 fOutput->Add(hnCandEv);
696
697 // Invariant mass related histograms
698 const Int_t nbinsmass = 200;
8577f7bc 699 const Int_t ptbinsD=100;
700 Float_t ptmin=0.,ptmax=50.;
701 TH2F *hInvMass = new TH2F("hInvMassptD", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, ptbinsD, ptmin, ptmax);
c683985a 702 hInvMass->SetStats(kTRUE);
703 hInvMass->GetXaxis()->SetTitle("mass (GeV/c)");
704 hInvMass->GetYaxis()->SetTitle("p_{T} (GeV/c)");
705 fOutput->Add(hInvMass);
8577f7bc 706 if ((fCandidateType==kDstartoKpipi) || fUseMCInfo){
c683985a 707 TH1F* hnSBCandEv=new TH1F("hnSBCandEv", "Number of side bands candidates per event (after cuts);# cand/ev", 100, 0.,100.);
708 fOutput->Add(hnSBCandEv);
8577f7bc 709 }
710 if (fCandidateType==kDstartoKpipi) {
c683985a 711
712 TH1F* hPtPion = new TH1F("hPtPion", "Primary pions candidates pt", 500, 0., 10.);
713 hPtPion->SetStats(kTRUE);
714 hPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");
715 hPtPion->GetYaxis()->SetTitle("entries");
716 fOutput->Add(hPtPion);
717 }
718
8577f7bc 719 const Int_t nbinsalpha=200;
720 Float_t minalpha=-TMath::Pi(), maxalpha=TMath::Pi();
721 const Int_t nbinsdeltaR= 200;
722 Float_t mindeltaR = 0., maxdeltaR = 10.;
723 if(fUseMCInfo){
724 if (fCandidateType==kDstartoKpipi){
725 TH2F* halphaDDS =new TH2F("halphaDDS","Angle D^{*}-D^{0} (Signal);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
726 TH2F* halphaDpisS=new TH2F("halphaDpisS","Angle D^{*}-#pi_{soft} (Signal);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
727 TH2F* halphaDpiS =new TH2F("halphaDpiS","Angle D^{*}-#pi (Signal);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
728 TH2F* halphaDKS =new TH2F("halphaDKS","Angle D^{*}-K (Signal);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
729
730 TH2F* halphaDDB =new TH2F("halphaDDB","Angle D^{*}-D^{0} (Background);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
731 TH2F* halphaDpisB=new TH2F("halphaDpisB","Angle D^{*}-#pi_{soft} (Background);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
732 TH2F* halphaDpiB =new TH2F("halphaDpiB","Angle D^{*}-#pi (Background);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
733 TH2F* halphaDKB =new TH2F("halphaDKB","Angle D^{*}-K (Background);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
734
735 TH2F* hdeltaRDDS =new TH2F("hdeltaRDDS","Angle D^{*}-D^{0} (Signal);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
736 TH2F* hdeltaRDpisS=new TH2F("hdeltaRDpisS","Angle D^{*}-#pi_{soft} (Signal);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
737 TH2F* hdeltaRDpiS =new TH2F("hdeltaRDpiS","Angle D^{*}-#pi (Signal);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
738 TH2F* hdeltaRDKS =new TH2F("hdeltaRDKS","Angle D^{*}-K (Signal);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
739
740 TH2F* hdeltaRDDB =new TH2F("hdeltaRDDB","Angle D^{*}-D^{0} (Background);#varphi (D^{*}) - #varphi (D0);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
741 TH2F* hdeltaRDpisB=new TH2F("hdeltaRDpisB","Angle D^{*}-#pi_{soft} (Background);#varphi (D^{*}) - #varphi (#pi_{soft});p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
742 TH2F* hdeltaRDpiB =new TH2F("hdeltaRDpiB","Angle D^{*}-#pi (Background);#varphi (D^{*}) - #varphi (#pi);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
743 TH2F* hdeltaRDKB =new TH2F("hdeltaRDKB","Angle D^{*}-K (Background);#varphi (D^{*}) - #varphi (K);p_{T}^{D*}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
744
745 fOutput->Add(halphaDDS);
746 fOutput->Add(halphaDpisS);
747 fOutput->Add(halphaDpiS);
748 fOutput->Add(halphaDKS);
749 fOutput->Add(halphaDDB);
750 fOutput->Add(halphaDpisB);
751 fOutput->Add(halphaDpiB);
752 fOutput->Add(halphaDKB);
753
754 fOutput->Add(hdeltaRDDS);
755 fOutput->Add(hdeltaRDpisS);
756 fOutput->Add(hdeltaRDpiS);
757 fOutput->Add(hdeltaRDKS);
758 fOutput->Add(hdeltaRDDB);
759 fOutput->Add(hdeltaRDpisB);
760 fOutput->Add(hdeltaRDpiB);
761 fOutput->Add(hdeltaRDKB);
762 }
763
764 if (fCandidateType==kD0toKpi){
765
766 TH2F* halphaDpiS=new TH2F("halphaDpiS","Angle D^{0}-#pi (Signal);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
767 TH2F* halphaDKS =new TH2F("halphaDKS","Angle D^{0}-K (Signal);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
768 TH2F* halphaDpiR=new TH2F("halphaDpiR","Angle D^{0}-#pi (Reflections);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
769 TH2F* halphaDKR =new TH2F("halphaDKR","Angle D^{0}-K (Reflections);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
770
771 TH2F* halphaDpiB=new TH2F("halphaDpiB","Angle D^{0}-#pi (Background);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
772 TH2F* halphaDKB =new TH2F("halphaDKB","Angle D^{0}-K (Background);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsalpha, minalpha, maxalpha, ptbinsD, ptmin, ptmax);
773
774
775 TH2F* hdeltaRDpiS=new TH2F("hdeltaRDpiS","Angle D^{0}-#pi (Signal);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
776 TH2F* hdeltaRDKS =new TH2F("hdeltaRDKS","Angle D^{0}-K (Signal);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
777 TH2F* hdeltaRDpiR=new TH2F("hdeltaRDpiR","Angle D^{0}-#pi (Reflections);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
778 TH2F* hdeltaRDKR =new TH2F("hdeltaRDKR","Angle D^{0}-K (Reflections);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
779
780 TH2F* hdeltaRDpiB=new TH2F("hdeltaRDpiB","Angle D^{0}-#pi (Background);#varphi (D^{0}) - #varphi (#pi);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
781 TH2F* hdeltaRDKB =new TH2F("hdeltaRDKB","Angle D^{0}-K (Background);#varphi (D^{0}) - #varphi (K);p_{T}^{D0}",nbinsdeltaR, mindeltaR, maxdeltaR, ptbinsD, ptmin, ptmax);
782
783 fOutput->Add(halphaDpiS);
784 fOutput->Add(halphaDKS);
785 fOutput->Add(halphaDpiR);
786 fOutput->Add(halphaDKR);
787 fOutput->Add(halphaDpiB);
788 fOutput->Add(halphaDKB);
789
790 fOutput->Add(hdeltaRDpiS);
791 fOutput->Add(hdeltaRDKS);
792 fOutput->Add(hdeltaRDpiR);
793 fOutput->Add(hdeltaRDKR);
794 fOutput->Add(hdeltaRDpiB);
795 fOutput->Add(hdeltaRDKB);
796
797 }
798
799 }
c683985a 800 return kTRUE;
801}
8577f7bc 802
803//_______________________________________________________________________________
804
805Float_t AliAnalysisTaskSEDmesonsFilterCJ::DeltaR(AliVParticle *p1, AliVParticle *p2) const {
806 //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
807
808 if(!p1 || !p2) return -1;
809 Double_t phi1=p1->Phi(),eta1=p1->Eta();
810 Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
811
812 Double_t dPhi=phi1-phi2;
813 if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
814 if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
815
816 Double_t dEta=eta1-eta2;
817 Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
818 return deltaR;
819
820}