Initial commit of the tasks for the flavor jet analysis; targeted at D-meson+jet...
[u/mrichter/AliRoot.git] / PWGJE / FlavourJetTasks / AliAnalysisTaskSEDmesonsFilterCJ.cxx
1 // $Id$
2
3 /**************************************************************************
4  * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
5  *                                                                        *
6  * Author: The ALICE Off-line Project.                                    *
7  * Contributors are mentioned in the code where appropriate.              *
8  *                                                                        *
9  * Permission to use, copy, modify and distribute this software and its   *
10  * documentation strictly for non-commercial purposes is hereby granted   *
11  * without fee, provided that the above copyright notice appears in all   *
12  * copies and that both the copyright notice and this permission notice   *
13  * appear in the supporting documentation. The authors make no claims     *
14  * about the suitability of this software for any purpose. It is          *
15  * provided "as is" without express or implied warranty.                  *
16  **************************************************************************/
17
18 /////////////////////////////////////////////////////////////
19 //
20 // AliAnalysisTaskSE for Dmesons - jet correlations analysis
21 //
22 // Author: Xiaoming Zhang, xmzhang@lbl.gov
23 ///////////////////////////////////////////////////////////////
24
25 #include <TMath.h>
26 #include <TH1D.h>
27 #include <TH2D.h>
28 #include <TList.h>
29 #include <TClonesArray.h>
30 #include <TDatabasePDG.h>
31
32 #include "AliAnalysisManager.h"
33 #include "AliAODTrack.h"
34 #include "AliAODEvent.h"
35 #include "AliAODHandler.h"
36 #include "AliAODExtension.h"
37 #include "AliAODRecoDecayHF.h"
38 #include "AliAODRecoDecayHF2Prong.h"
39 #include "AliAODRecoCascadeHF.h"
40 #include "AliRDHFCuts.h"
41 #include "AliRDHFCutsD0toKpi.h"
42 #include "AliRDHFCutsDStartoKpipi.h"
43 #include "AliAnalysisTaskSEDmesonsFilterCJ.h"
44
45 ClassImp(AliAnalysisTaskSEDmesonsFilterCJ)
46
47 //_____________________________________________________________________________
48 AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ() :
49 AliAnalysisTaskSE(),
50 fEventAOD(0),
51 fCutDzero(0),
52 fCutDstar(0),
53 fDzeroClArr(0),
54 fDstarClArr(0),
55 fUsedDzeros(0),
56 fUsedD0bars(0),
57 fUsedDstars(0),
58 fIsNotExecOnce(kTRUE),
59 fListCutsDmesons(0),
60 fListDzeroHistos(0),
61 fListDstarHistos(0)
62 {
63 //
64 // Default constructor
65 //
66 }
67
68 //_____________________________________________________________________________
69 AliAnalysisTaskSEDmesonsFilterCJ::AliAnalysisTaskSEDmesonsFilterCJ(const char *name) :
70 AliAnalysisTaskSE(name),
71 fEventAOD(0),
72 fCutDzero(0),
73 fCutDstar(0),
74 fDzeroClArr(0),
75 fDstarClArr(0),
76 fUsedDzeros(0),
77 fUsedD0bars(0),
78 fUsedDstars(0),
79 fIsNotExecOnce(kTRUE),
80 fListCutsDmesons(0),
81 fListDzeroHistos(0),
82 fListDstarHistos(0)
83 {
84 //
85 // Constructor
86 //
87
88   DefineOutput(1, TList::Class());
89   DefineOutput(2, TList::Class());
90   DefineOutput(3, TList::Class());
91 }
92
93 //_____________________________________________________________________________
94 AliAnalysisTaskSEDmesonsFilterCJ::~AliAnalysisTaskSEDmesonsFilterCJ()
95 {
96 //
97 // Default destructor
98 //
99
100   if (fEventAOD)   { delete fEventAOD;   fEventAOD  =0; }
101
102   if (fCutDzero)   { delete fCutDzero;   fCutDzero  =0; }
103   if (fCutDstar)   { delete fCutDstar;   fCutDstar  =0; }
104
105   if (fDzeroClArr) { delete fDzeroClArr; fDzeroClArr=0; }
106   if (fDstarClArr) { delete fDstarClArr; fDstarClArr=0; }
107
108   if (fUsedDzeros) { delete fUsedDzeros; fUsedDzeros=0; }
109   if (fUsedD0bars) { delete fUsedD0bars; fUsedD0bars=0; }
110   if (fUsedDstars) { delete fUsedDstars; fUsedDstars=0; }
111
112   if (fListCutsDmesons) { delete fListCutsDmesons; fListCutsDmesons=0; }
113   if (fListDzeroHistos) { delete fListDzeroHistos; fListDzeroHistos=0; }
114   if (fListDstarHistos) { delete fListDstarHistos; fListDstarHistos=0; }
115 }
116
117 //_____________________________________________________________________________
118 void AliAnalysisTaskSEDmesonsFilterCJ::Init()
119 {
120 //
121 // AliAnalysisTaskSEDmesonsFilterCJ::Init()
122 //
123
124   if (fListCutsDmesons) return;
125   fListCutsDmesons = new TList(); fListCutsDmesons->SetOwner();
126
127   if (fCutDzero) {
128     AliRDHFCutsD0toKpi *cutDzero = new AliRDHFCutsD0toKpi(*fCutDzero);
129     cutDzero->SetName("AnalysisCutsDzero"); 
130     fListCutsDmesons->Add(cutDzero);
131   }
132
133   if (fCutDstar) {
134     AliRDHFCutsDStartoKpipi *cutDstar = new AliRDHFCutsDStartoKpipi(*fCutDstar);
135     cutDstar->SetName("AnalysisCutsDstar"); 
136     fListCutsDmesons->Add(cutDstar);
137   }
138
139   return;
140 }
141
142 //_____________________________________________________________________________
143 void AliAnalysisTaskSEDmesonsFilterCJ::UserCreateOutputObjects()
144 {
145 //
146 // AliAnalysisTaskSEDmesonsFilterCJ::UserCreateOutputObjects
147 //
148
149   if (fCutDzero) { fUsedDzeros = new TClonesArray("AliAODRecoDecayHF2Prong",0); fUsedDzeros->SetName("AnaUsedDzero");}
150   if (fCutDzero) { fUsedD0bars = new TClonesArray("AliAODRecoDecayHF2Prong",0); fUsedD0bars->SetName("AnaUsedD0bar");}
151   if (fCutDstar) { fUsedDstars = new TClonesArray("AliAODRecoCascadeHF",    0); fUsedDstars->SetName("AnaUsedDstar"); }
152
153   if (!fListDzeroHistos) fListDzeroHistos = new TList(); fListDzeroHistos->SetOwner();
154   if (!fListDstarHistos) fListDstarHistos = new TList(); fListDstarHistos->SetOwner();
155
156   if (fCutDzero) CreateDzeroHistograms();
157   if (fCutDstar) CreateDstarHistograms();
158
159   PostData(1, fListCutsDmesons);
160   PostData(2, fListDzeroHistos);
161   PostData(3, fListDstarHistos);
162   return;
163 }
164
165 //_____________________________________________________________________________
166 void AliAnalysisTaskSEDmesonsFilterCJ::UserExec(Option_t *)
167 {
168 //
169 // AliAnalysisTaskSEDmesonsFilterCJ::UserExec
170 //
171
172   if (fIsNotExecOnce) {
173     fIsNotExecOnce=ExecOnce();
174     if (fIsNotExecOnce) return;
175   }
176
177   if (ExecEach()) ExecAnas();
178
179   return;
180 }
181
182 //_____________________________________________________________________________
183 void AliAnalysisTaskSEDmesonsFilterCJ::Terminate(Option_t *)
184 {
185 //
186 // AliAnalysisTaskSEDmesonsFilterCJ::Terminate
187 //
188
189   return;
190 }
191
192 //_____________________________________________________________________________
193 void AliAnalysisTaskSEDmesonsFilterCJ::NotifyRun()
194 {
195 //
196 // AliAnalysisTaskSEDmesonsFilterCJ::NotifyRun
197 //
198
199   return;
200 }
201
202 //_____________________________________________________________________________
203 void AliAnalysisTaskSEDmesonsFilterCJ::ExecAnas()
204 {
205 //
206 // AliAnalysisTaskSEDmesonsFilterCJ::ExecAnas
207 //
208
209   if (fUsedDzeros) fUsedDzeros->Delete();
210   if (fUsedD0bars) fUsedD0bars->Delete();
211   if (fUsedDstars) fUsedDstars->Delete();
212
213   if (fDzeroClArr) ExecAnasDzero();
214   if (fDstarClArr) ExecAnasDstar();
215   return;
216 }
217
218 //_____________________________________________________________________________
219 void AliAnalysisTaskSEDmesonsFilterCJ::ExecAnasDzero()
220 {
221 //
222 // AliAnalysisTaskSEDmesonsFilterCJ::ExecAnasDzero
223 //
224
225   if (!fCutDzero->IsEventSelected(fEventAOD)) return;
226   const Int_t nCands = fDzeroClArr->GetEntriesFast(); if (nCands==0) return;
227
228   ((TH1D*)fListDzeroHistos->FindObject("hDzero_candi"))->Fill(nCands);
229
230   Int_t countN = 0;
231   AliAODRecoDecayHF2Prong *candi = 0;
232   for (Int_t iCand=0; iCand<nCands; iCand++) {
233     candi = dynamic_cast<AliAODRecoDecayHF2Prong*>(fDzeroClArr->At(iCand)); if (!candi) continue;
234
235     Double_t dPt = candi->Pt();
236     if (!fCutDzero->IsInFiducialAcceptance(dPt,candi->YD0())) continue;
237     Int_t mask = fCutDzero->IsSelected(candi,AliRDHFCuts::kAll,fEventAOD); if (mask==0) continue;
238
239     Double_t dMass = 0.; 
240     if ((mask==1) || (mask==3)) dMass = candi->InvMassD0();
241     if  (mask==2)               dMass = candi->InvMassD0bar();
242     ((TH2D*)fListDzeroHistos->FindObject("hDzero_Mass_Pt"))->Fill(dMass,dPt);
243
244     if ((mask==1) || (mask==3)) new((*fUsedDzeros)[countN++]) AliAODRecoDecayHF2Prong(*candi);
245     if  (mask==2)               new((*fUsedD0bars)[countN++]) AliAODRecoDecayHF2Prong(*candi);
246     candi = 0;
247   }
248
249   ((TH1D*)fListDzeroHistos->FindObject("hDzero_seled"))->Fill(countN);
250   return;
251 }
252
253 //_____________________________________________________________________________
254 void AliAnalysisTaskSEDmesonsFilterCJ::ExecAnasDstar()
255 {
256 //
257 //  AliAnalysisTaskSEDmesonsFilterCJ::ExecAnasDstar
258 //
259
260   if (!fCutDstar->IsEventSelected(fEventAOD)) return;
261   const Int_t nCands = fDstarClArr->GetEntriesFast(); if (nCands==0) return;
262
263   Double_t dMassDzero = TDatabasePDG::Instance()->GetParticle(421)->Mass();
264   Double_t dMassDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
265   Double_t dMassDelta = dMassDstar - dMassDzero;
266   Double_t dSigmaD0Pt[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 };
267
268   ((TH1D*)fListDstarHistos->FindObject("hDstar_candi"))->Fill(nCands);
269
270   Int_t countN = 0;
271   AliAODRecoCascadeHF *candi = 0;
272   for (Int_t iCand=0; iCand<nCands; iCand++) {
273     candi = dynamic_cast<AliAODRecoCascadeHF*>(fDstarClArr->At(iCand)); if (!candi) continue;
274
275     Double_t dPt = candi->Pt();
276     if (!fCutDstar->IsInFiducialAcceptance(dPt,candi->YDstar()))   continue;
277     if (!fCutDstar->IsSelected(candi,AliRDHFCuts::kAll)) continue;
278
279     Int_t iBin = fCutDzero->PtBin(dPt);
280     if ((iBin<0) || (iBin>=fCutDzero->GetNPtBins())) { AliWarning(Form("pT(D*)=%f out of bounds!!!",dPt)); continue; }
281
282     Double_t dMassResD0 = TMath::Abs(candi->InvMassD0() - dMassDzero); if (dMassResD0>3.*dSigmaD0Pt[iBin]) continue;
283
284     Double_t deltaM = candi->DeltaInvMass();
285     ((TH2D*)fListDstarHistos->FindObject("hDstar_deltaMass_Pt"))->Fill(deltaM,dPt);
286
287     if (TMath::Abs(deltaM-dMassDelta)) ((TH1D*)fListDstarHistos->FindObject("hDstar_SoftPiPt"))->Fill(((AliAODTrack*)candi->GetBachelor())->Pt());
288
289     new((*fUsedDstars)[countN++]) AliAODRecoCascadeHF(*candi);
290     candi = 0;
291   }
292
293   ((TH1D*)fListDstarHistos->FindObject("hDstar_seled"))->Fill(countN);
294   return;
295 }
296
297 //_____________________________________________________________________________
298 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::ExecOnce()
299 {
300 //
301 // AliAnalysisTaskSEDmesonsFilterCJ::ExecOnce
302 //
303
304   if (!InputEvent())  return kTRUE;
305   fEventAOD = dynamic_cast<AliAODEvent*>(InputEvent());
306
307   if (fEventAOD) {
308     if (fCutDzero) fDzeroClArr = dynamic_cast<TClonesArray*>(fEventAOD->FindListObject("D0toKpi"));
309     if (fCutDstar) fDstarClArr = dynamic_cast<TClonesArray*>(fEventAOD->FindListObject("Dstar"));
310   } else if (AODEvent() && IsStandardAOD()) {
311     fEventAOD = dynamic_cast<AliAODEvent*>(AODEvent());
312     AliAODHandler* aodH = dynamic_cast<AliAODHandler*>((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
313
314     if(aodH->GetExtensions()) {
315       AliAODExtension    *aodExt = dynamic_cast<AliAODExtension*>(aodH->GetExtensions()->FindObject("AliAOD.VertexingHF.root"));
316       if (fCutDzero) fDzeroClArr = dynamic_cast<TClonesArray*>(aodExt->GetAOD()->FindListObject("D0toKpi"));
317       if (fCutDstar) fDstarClArr = dynamic_cast<TClonesArray*>(aodExt->GetAOD()->FindListObject("Dstar"));
318     }
319   }
320
321   if (!fEventAOD)                { AliError("No AOD Event!!!");    return kTRUE; }
322   if (!fDzeroClArr && fCutDzero) { AliError("No Dzero Branch!!!"); return kTRUE; }
323   if (!fDstarClArr && fCutDstar) { AliError("No Dstar Branch!!!"); return kTRUE; }
324
325   if (fDzeroClArr && fUsedDzeros) {
326     fUsedDzeros->Delete();
327     if (!(InputEvent()->FindListObject("AnaUsedDzero"))) InputEvent()->AddObject(fUsedDzeros);
328   }
329
330   if (fDzeroClArr && fUsedD0bars) {
331     fUsedD0bars->Delete();
332     if (!(InputEvent()->FindListObject("AnaUsedD0bar"))) InputEvent()->AddObject(fUsedD0bars);
333   }
334
335   if (fDstarClArr && fUsedDstars) {
336     fUsedDstars->Delete();
337     if (!(InputEvent()->FindListObject("AnaUsedDstar"))) InputEvent()->AddObject(fUsedDstars);
338   }
339
340   return kFALSE;
341 }
342
343 //_____________________________________________________________________________
344 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::ExecEach()
345 {
346 //
347 // AliAnalysisTaskSEDmesonsFilterCJ::ExecEach
348 //
349
350   if (!fEventAOD) {
351     AliWarning("No Input Event, Skip Event!!!");
352     return kFALSE;
353   }
354
355   if ((!fDzeroClArr) && (!fDstarClArr)) {
356     AliWarning("No D meson Branch, Skip Event!!!");
357     return kFALSE;
358   }
359
360   if ((!fEventAOD->GetPrimaryVertex()) || (TMath::Abs(fEventAOD->GetMagneticField())<1e-3)) return kFALSE;
361
362   Bool_t bIsDmeson = kTRUE;
363   if (fDzeroClArr) bIsDmeson = ((fDzeroClArr->GetEntriesFast()==0) && bIsDmeson);
364   if (fDstarClArr) bIsDmeson = ((fDstarClArr->GetEntriesFast()==0) && bIsDmeson);
365   if (bIsDmeson)   return kFALSE;
366
367   return kTRUE;
368 }
369
370 //_____________________________________________________________________________
371 void AliAnalysisTaskSEDmesonsFilterCJ::CreateDzeroHistograms()
372 {
373 //
374 // AliAnalysisTaskSEDmesonsFilterCJ::CreateDzeroHistograms
375 //
376
377   if (!fListDzeroHistos) return;
378   Bool_t bStatusTmpH = TH1::AddDirectoryStatus();
379   TH1::AddDirectory(kFALSE);
380
381   TH1D *h1D = 0;
382   h1D = new TH1D("hDzero_candi", "", 100, -0.5, 99.5); h1D->Sumw2(); fListDzeroHistos->Add(h1D); h1D=0;
383   h1D = new TH1D("hDzero_seled", "", 100, -0.5, 99.5); h1D->Sumw2(); fListDzeroHistos->Add(h1D); h1D=0;
384
385   TH2D *h2D = 0;
386   Double_t dMass = TDatabasePDG::Instance()->GetParticle(421)->Mass();
387   h2D = new TH2D("hDzero_Mass_Pt", "", 200, dMass-0.15, dMass+0.15, 100, 0., 50.);
388   h2D->Sumw2(); fListDzeroHistos->Add(h2D); h2D=0;
389
390   TH1::AddDirectory(bStatusTmpH);
391   return;
392 }
393
394 //_____________________________________________________________________________
395 void AliAnalysisTaskSEDmesonsFilterCJ::CreateDstarHistograms()
396 {
397 //
398 // AliAnalysisTaskSEDmesonsFilterCJ::CreateDstarHistograms
399 //
400
401   if (!fListDstarHistos) return;
402   Bool_t bStatusTmpH = TH1::AddDirectoryStatus();
403   TH1::AddDirectory(kFALSE);
404
405   TH1D *h1D = 0;
406   h1D = new TH1D("hDstar_candi",    "", 100, -0.5, 99.5); h1D->Sumw2(); fListDstarHistos->Add(h1D); h1D=0;
407   h1D = new TH1D("hDstar_seled",    "", 100, -0.5, 99.5); h1D->Sumw2(); fListDstarHistos->Add(h1D); h1D=0;
408   h1D = new TH1D("hDstar_SoftPiPt", "", 500,  0.,  10.0); h1D->Sumw2(); fListDstarHistos->Add(h1D); h1D=0;
409
410   TH2D *h2D = 0;
411   Double_t dMassDzero = TDatabasePDG::Instance()->GetParticle(421)->Mass();
412   Double_t dMassDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
413   Double_t dMassDelta = dMassDstar - dMassDzero;
414   h2D = new TH2D("hDstar_deltaMass_Pt", "", 200, dMassDelta-0.015, dMassDelta+0.015, 100, 0., 50.);
415   h2D->Sumw2(); fListDstarHistos->Add(h2D); h2D=0;
416
417   TH1::AddDirectory(bStatusTmpH);
418   return;
419 }
420
421 //_____________________________________________________________________________
422 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::SetCutDzero(AliRDHFCutsD0toKpi *cut)
423 {
424 //
425 // AliAnalysisTaskSEDmesonsFilterCJ::SetCutDzero
426 //
427
428   if (!cut) return kTRUE;
429
430   if (fCutDzero) { delete fCutDzero; fCutDzero=0; }
431   fCutDzero = new AliRDHFCutsD0toKpi(*cut);
432   if (!fCutDzero) { AliError("No Dzero Cuts!!!"); return kTRUE; }
433
434   return kFALSE;
435 }
436
437 //_____________________________________________________________________________
438 Bool_t AliAnalysisTaskSEDmesonsFilterCJ::SetCutDstar(AliRDHFCutsDStartoKpipi *cut)
439 {
440 //
441 // AliAnalysisTaskSEDmesonsFilterCJ::SetCutDstar
442 //
443
444   if (!cut) return kTRUE;
445
446   if (fCutDstar) { delete fCutDstar; fCutDstar=0; }
447   fCutDstar = new AliRDHFCutsDStartoKpipi(*cut);
448   if (!fCutDstar)  { AliError("No Dstar Cuts!!!"); return kTRUE; }
449
450   return kFALSE;
451 }