]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/charmFlow/AliAnalysisTaskSEHFv2.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
adjusted to use proper cluster timing for EMCAL
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / charmFlow / AliAnalysisTaskSEHFv2.cxx
CommitLineData
a8f6c03f 1/**************************************************************************
2 * Copyright(c) 1998-2010, 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// AliAnalysisTaskSEHFv2 gives the needed tools for the D
19// mesons v2 analysis with event plane method
20// Authors: Chiara Bianchin, cbianchi@pd.infn.it,
21// Robert Grajcarek, grajcarek@physi.uni-heidelberg.de
22// Giacomo Ortona, ortona@to.infn.it
23// Carlos Eugenio Perez Lara, carlos.eugenio.perez.lara@cern.ch
24//
25/////////////////////////////////////////////////////////////
26
f829c4f7 27/* $Id$ */
28
a8f6c03f 29#include <Riostream.h>
30#include <TClonesArray.h>
31#include <TCanvas.h>
32#include <TList.h>
33#include <TFile.h>
34#include <TString.h>
35#include <TH1F.h>
36#include <TH2F.h>
37#include <TGraphErrors.h>
38#include <TGraph.h>
39#include <TDatabasePDG.h>
40#include <TRandom3.h>
41#include <TVector2.h>
42#include <TArrayF.h>
43
44#include <AliLog.h>
45#include <AliAnalysisDataSlot.h>
46#include <AliAnalysisDataContainer.h>
47#include "AliAnalysisManager.h"
48#include "AliAODHandler.h"
49#include "AliAODEvent.h"
50#include "AliAODVertex.h"
51#include "AliAODTrack.h"
52#include "AliAODMCHeader.h"
53#include "AliAODMCParticle.h"
54#include "AliAODRecoDecayHF3Prong.h"
55#include "AliAODRecoDecayHF.h"
56#include "AliAODRecoDecayHF2Prong.h"
57#include "AliAODRecoDecayHF4Prong.h"
58#include "AliAODRecoCascadeHF.h"
59#include "AliAODHFUtil.h"
60
61#include "AliAnalysisTaskSE.h"
62#include "AliRDHFCutsDplustoKpipi.h"
63#include "AliRDHFCutsD0toKpipipi.h"
64#include "AliRDHFCutsDstoKKpi.h"
65#include "AliRDHFCutsDStartoKpipi.h"
66#include "AliRDHFCutsD0toKpi.h"
67#include "AliRDHFCutsLctopKpi.h"
68
69#include "AliHFMassFitter.h"
70#include "AliEventplane.h"
71#include "AliFlowTrack.h"
72#include "AliFlowVector.h"
73#include "AliFlowTrackCuts.h"
74#include "AliFlowEvent.h"
75
76#include "AliAnalysisTaskSEHFv2.h"
77
78ClassImp(AliAnalysisTaskSEHFv2)
79
80
81//________________________________________________________________________
82AliAnalysisTaskSEHFv2::AliAnalysisTaskSEHFv2():
83AliAnalysisTaskSE(),
84 fhEventsInfo(0),
85 fOutput(0),
86 fRDCuts(0),
87 fParHist(0),
88 fLowmasslimit(1.765),
89 fUpmasslimit(1.965),
90 fNPtBins(1),
91 fNPhiBinLims(2),
92 fPhiBins(0),
93 fCentLowLimit(0),
94 fCentUpLimit(100),
95 fNMassBins(200),
96 fReadMC(kFALSE),
97 fDecChannel(0),
98 fUseV0EP(kFALSE),
99 fV0EPorder(2)
100{
101 // Default constructor
102}
103
104//________________________________________________________________________
105AliAnalysisTaskSEHFv2::AliAnalysisTaskSEHFv2(const char *name,AliRDHFCuts *rdCuts,Int_t decaychannel, Int_t nlimsphibin, Float_t *phibinlimits,TH2D **histPar):
106 AliAnalysisTaskSE(name),
107 fhEventsInfo(0),
108 fOutput(0),
109 fRDCuts(rdCuts),
110 fParHist(0),
111 fLowmasslimit(0),
112 fUpmasslimit(0),
113 fNPtBins(1),
114 fNPhiBinLims(2),
115 fPhiBins(0),
116 fCentLowLimit(0),
117 fCentUpLimit(100),
118 fNMassBins(200),
119 fReadMC(kFALSE),
120 fDecChannel(decaychannel),
121 fUseV0EP(kFALSE),
122 fV0EPorder(2)
123{
124
125 Int_t pdg=421;
126 switch(fDecChannel){
127 case 0:
128 pdg=411;
129 break;
130 case 1:
131 pdg=421;
132 break;
133 case 2:
134 pdg=413;
135 break;
136 case 3:
137 pdg=431;
138 break;
139 case 4:
140 pdg=421;
141 break;
142 case 5:
143 pdg=4122;
144 break;
145 }
146 for(Int_t i=0;i<6;i++)fHistvzero[i]=(TH2D*)histPar[i]->Clone();
147 for(Int_t i=0;i<6;i++)if(!fHistvzero[i])printf("No VZERO histograms!\n");
148 if(pdg==413) SetMassLimits((Float_t)0.135,(Float_t)0.165);
149 else SetMassLimits((Float_t)0.2,pdg); //check range
150 fNPtBins=fRDCuts->GetNPtBins();
151 if(nlimsphibin>2) fNPhiBinLims=nlimsphibin;
152 else AliInfo("At least 2 limits in Delta phi needed");
153
154 fPhiBins=new Float_t[fNPhiBinLims];
155 for(Int_t i=0;i<fNPhiBinLims;i++) fPhiBins[i]=phibinlimits[i];
156
157 if(fDebug>1)fRDCuts->PrintAll();
158 // Output slot #1 writes into a TH1F container
159 DefineOutput(1,TH1F::Class()); //Info on the number of events etc.
160 // Output slot #2 writes into a TList container
161 DefineOutput(2,TList::Class()); //Main output
162 // Output slot #3 writes into a AliRDHFCuts container (cuts)
163 switch(fDecChannel){
164 case 0:
165 DefineOutput(3,AliRDHFCutsDplustoKpipi::Class()); //Cut object for Dplus
166 break;
167 case 1:
168 DefineOutput(3,AliRDHFCutsD0toKpi::Class()); //Cut object for D0
169 break;
170 case 2:
171 DefineOutput(3,AliRDHFCutsDStartoKpipi::Class()); //Cut object for D*
172 break;
173 }
174 //DefineOutput(4,AliFlowEventSimple::Class());
175 DefineOutput(4,TList::Class());
176}
177
178//________________________________________________________________________
179AliAnalysisTaskSEHFv2::~AliAnalysisTaskSEHFv2()
180{
181 // Destructor
182 if (fOutput) {
183 delete fOutput;
184 fOutput = 0;
185 }
186
187 if(fhEventsInfo){
188 delete fhEventsInfo;
189 fhEventsInfo=0;
190 }
191
192 if(fRDCuts) {
193 delete fRDCuts;
194 fRDCuts= 0;
195 }
196
197 if(fParHist) {
198 delete fParHist;
199 fParHist= 0;
200 }
201 for(Int_t i=0;i<6;i++){
202 if(fHistvzero[i]) {
203 delete fHistvzero[i];
204 fHistvzero[i]=0x0;
205 }
206 }
207
208}
209
210//_________________________________________________________________
211void AliAnalysisTaskSEHFv2::SetMassLimits(Float_t range, Int_t pdg){
212 Float_t mass=0;
213 Int_t abspdg=TMath::Abs(pdg);
214 mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
215 fUpmasslimit = mass+range;
216 fLowmasslimit = mass-range;
217}
218//_________________________________________________________________
219void AliAnalysisTaskSEHFv2::SetMassLimits(Float_t lowlimit, Float_t uplimit){
220 if(uplimit>lowlimit)
221 {
222 fUpmasslimit = uplimit;
223 fLowmasslimit = lowlimit;
224 }
225}
226
227
228//________________________________________________________________________
229void AliAnalysisTaskSEHFv2::LocalInit()
230{
231 // Initialization
232
233 if(fDebug > 1) printf("AnalysisTaskSEHFv2::Init() \n");
234
235 fRDCuts->SetMinCentrality(10);
236 fRDCuts->SetMaxCentrality(80);
237
238 switch(fDecChannel){
239 case 0:
240 {
241 AliRDHFCutsDplustoKpipi* copycut=new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fRDCuts)));
242 // Post the data
243 PostData(3,copycut);
244 }
245 break;
246 case 1:
247 {
248 AliRDHFCutsD0toKpi* copycut=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fRDCuts)));
249 // Post the data
250 PostData(3,copycut);
251 }
252 break;
253 case 2:
254 {
255 AliRDHFCutsDStartoKpipi* copycut=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fRDCuts)));
256 // Post the data
257 PostData(3,copycut);
258 }
259 break;
260 default:
261 return;
262 }
263 return;
264}
265//________________________________________________________________________
266void AliAnalysisTaskSEHFv2::UserCreateOutputObjects()
267{
268 // Create the output container
269
270 if(fDebug > 1) printf("AnalysisTaskSEHFv2::UserCreateOutputObjects() \n");
271
272 fhEventsInfo=new TH1F(GetOutputSlot(1)->GetContainer()->GetName(), "Number of AODs scanned",7,-0.5,6.5);
273 fhEventsInfo->GetXaxis()->SetBinLabel(1,"nEventsAnal");
274 fhEventsInfo->GetXaxis()->SetBinLabel(2,"nEvSelected");
275 fhEventsInfo->GetXaxis()->SetBinLabel(3,"nCandidatesSelected");
276 fhEventsInfo->GetXaxis()->SetBinLabel(4,"out of pt bounds");
277 fhEventsInfo->GetXaxis()->SetBinLabel(5,"Pile-up Rej");
278 fhEventsInfo->GetXaxis()->SetBinLabel(6,"N. of 0SMH");
279 fhEventsInfo->GetXaxis()->SetBinLabel(7,Form("Ev Sel in Centr %.0f-%.0f%s",fRDCuts->GetMinCentrality(),fRDCuts->GetMaxCentrality(),"%"));
280 fhEventsInfo->GetXaxis()->SetNdivisions(1,kFALSE);
281
282
283 // Several histograms are more conveniently managed in a TList
284 fOutput = new TList();
285 fOutput->SetOwner();
286 fOutput->SetName("MainOutput");
287
288 for(Int_t icentr=10;icentr<=80;icentr=icentr+5){
289 TString centrname;centrname.Form("centr%d_%d",icentr-5,icentr);
290 for(Int_t i=0;i<fNPtBins;i++){
291
292 TString hname;
293 TString title;
294 hname.Form("hPhi_pt%d",i);hname.Append(centrname.Data());
295 title.Form("Phi distribution (Pt bin %d %s);#phi;Entries",i,centrname.Data());
296
297 TH1F* hPhi=new TH1F(hname.Data(),title.Data(),96,0.,2*TMath::Pi());
298 hPhi->Sumw2();
299 fOutput->Add(hPhi);
300
301 for(Int_t j=0;j<fNPhiBinLims-1;j++){
302
303 hname.Form("hMass_pt%dphi%d",i,j);hname.Append(centrname.Data());
304 title.Form("Invariant mass (Pt bin %d, Phi bin %d, %s);M(GeV/c^{2});Entries",i,j,centrname.Data());
305
306 TH1F* hMass=new TH1F(hname.Data(),title.Data(),fNMassBins,fLowmasslimit,fUpmasslimit);
307 hMass->Sumw2();
308 fOutput->Add(hMass);
309
310
311 if(fReadMC){
312 hname.Form("hSgn_pt%dphi%d",i,j);hname.Append(centrname.Data());
313 title.Form("Invariant mass S (Pt bin %d, Phi bin %d, %s);M(GeV/c^{2});Entries",i,j,centrname.Data());
314 TH1F* hSgn=new TH1F(hname.Data(),title.Data(),fNMassBins,fLowmasslimit,fUpmasslimit);
315 hSgn->Sumw2();
316 fOutput->Add(hSgn);
317
318 hname.Form("hBkg_pt%dphi%d",i,j);hname.Append(centrname.Data());
319 title.Form("Invariant mass B (Pt bin %d, Phi bin %d, %s);M(GeV/c^{2});Entries",i,j,centrname.Data());
320 TH1F* hBkg=new TH1F(hname.Data(),title.Data(),fNMassBins,fLowmasslimit,fUpmasslimit);
321 hBkg->Sumw2();
322 fOutput->Add(hBkg);
323
324 if((fDecChannel != AliAnalysisTaskSEHFv2::kDplustoKpipi) && (fDecChannel != AliAnalysisTaskSEHFv2::kDstartoKpipi)){
325 hname.Form("hRfl_pt%dphi%d",i,j);hname.Append(centrname.Data());
326 title.Form("Invariant mass Reflections (Pt bin %d, Phi bin %d %s);M(GeV/c^{2});Entries",i,j,centrname.Data());
327 TH1F* hRfl=new TH1F(hname.Data(),title.Data(),fNMassBins,fLowmasslimit,fUpmasslimit);
328 hRfl->Sumw2();
329 fOutput->Add(hRfl);
330 }
331 }
332 }
333
334 TH2F* hMc2phi=new TH2F(Form("hMc2phi_pt%d%s",i,centrname.Data()),Form("Mass vs cos2#Delta#phi (p_{t} bin %d %s);cos2#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
335 fOutput->Add(hMc2phi);
336
337 if (fReadMC){
338 TH2F* hMc2phiS=new TH2F(Form("hMc2phiS_pt%d%s",i,centrname.Data()),Form("Mass vs cos2#Delta#phi (p_{t} bin %d %s);cos2#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
339 fOutput->Add(hMc2phiS);
340 TH2F * hMphiS=new TH2F(Form("hMphiS_pt%d%s",i,centrname.Data()),Form("Mass vs #Delta#phi (p_{t} bin %d %s);#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),96,0,2*TMath::Pi(),fNMassBins,fLowmasslimit,fUpmasslimit);
341 fOutput->Add(hMphiS);
342 TH2F* hMc2phiB=new TH2F(Form("hMc2phiB_pt%d%s",i,centrname.Data()),Form("Mass vs cos2#Delta#phi (p_{t} bin %d %s);cos2#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
343 fOutput->Add(hMc2phiB);
344 TH2F * hMphiB=new TH2F(Form("hMphiB_pt%d%s",i,centrname.Data()),Form("Mass vs #Delta#phi (p_{t} bin %d %s);#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),96,0,2*TMath::Pi(),fNMassBins,fLowmasslimit,fUpmasslimit);
345 fOutput->Add(hMphiB);
346 if((fDecChannel != AliAnalysisTaskSEHFv2::kDplustoKpipi) &&(fDecChannel != AliAnalysisTaskSEHFv2::kDstartoKpipi)){
347 TH2F* hMc2phiR=new TH2F(Form("hMc2phiR_pt%d%s",i,centrname.Data()),Form("Mass vs cos2#Delta#phi (p_{t} bin %d %s);cos2#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),100,-1.,1.,fNMassBins,fLowmasslimit,fUpmasslimit);
348 fOutput->Add(hMc2phiR);
349 TH2F* hMphiR=new TH2F(Form("hMphiR_pt%d%s",i,centrname.Data()),Form("Mass vs #Delta#phi (p_{t} bin %d %s);#Delta#phi;M (GeV/c^{2})",i,centrname.Data()),96,0,2*TMath::Pi(),fNMassBins,fLowmasslimit,fUpmasslimit);
350 fOutput->Add(hMphiR);
351 }
352 }
353 }
354
355 TH2F* hMphi=new TH2F(Form("hMphi%s",centrname.Data()),Form("Mass vs #Delta#phi %s;#Delta#phi;M (GeV/c^{2})",centrname.Data()),96,0,TMath::Pi(),fNMassBins,fLowmasslimit,fUpmasslimit);
356 fOutput->Add(hMphi);
357
358 TH1F* hEvPlane=new TH1F(Form("hEvPlane%s",centrname.Data()),Form("Event plane angle %s;#phi Ev Plane;Entries",centrname.Data()),200,0.,TMath::Pi());
359 fOutput->Add(hEvPlane);
360
361 //TH1F* hEvPlaneCheck=new TH1F(Form("hEvPlaneCheck%s",centrname.Data()),Form("Event plane angle - Event plane angle per candidate %s;(#phi(Ev Plane) - #phi(Ev Plane Candidate))/#phi(EvPlane);Entries",centrname.Data()),200,-0.2,0.2);
362 //fOutput->Add(hEvPlaneCheck);
363
364 TH1F* hEvPlaneCand=new TH1F(Form("hEvPlaneCand%s",centrname.Data()),Form("Event plane angle - Event plane angle per candidate %s;#phi(Ev Plane Candidate);Entries",centrname.Data()),200,-TMath::Pi(),TMath::Pi());
365 fOutput->Add(hEvPlaneCand);
366
367 TH1F* hEvPlaneReso=new TH1F(Form("hEvPlaneReso%s",centrname.Data()),Form("Event plane angle Resolution %s;cos2(#psi_{A}-#psi_{B});Entries",centrname.Data()),220,-1.1,1.1);
368 fOutput->Add(hEvPlaneReso);
369 }
370
371 TH1F* hPhiBins=new TH1F("hPhiBins","Bins in #Delta#phi used in this analysis;#phi bin;n jobs",fNPhiBinLims-1,fPhiBins);
372 fOutput->Add(hPhiBins);
373 for(Int_t k=0;k<fNPhiBinLims-1;k++)hPhiBins->SetBinContent(k+1,1);
374
375 PostData(1,fhEventsInfo);
376 PostData(2,fOutput);
377 fParHist = new TList();
378 fParHist->SetOwner();
379 fParHist->SetName("VZEROcorr");
380 for(Int_t i=0;i<6;i++){
381 fParHist->Add((TH2D*)fHistvzero[i]);
382 }
383 PostData(4,fParHist);
384 return;
385}
386
387//________________________________________________________________________
388void AliAnalysisTaskSEHFv2::UserExec(Option_t */*option*/)
389{
390 // Execute analysis for current event:
391 // heavy flavor candidates association to MC truth
392
393 AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
394 if(fDebug>2) printf("Analysing decay %d\n",fDecChannel);
395 // Post the data already here
396 PostData(1,fhEventsInfo);
397 PostData(2,fOutput);
398
399 TClonesArray *arrayProng =0;
400 Int_t absPdgMom=0;
401 if(!aod && AODEvent() && IsStandardAOD()) {
402 // In case there is an AOD handler writing a standard AOD, use the AOD
403 // event in memory rather than the input (ESD) event.
404 aod = dynamic_cast<AliAODEvent*> (AODEvent());
405 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
406 // have to taken from the AOD event hold by the AliAODExtension
407 AliAODHandler* aodHandler = (AliAODHandler*)
408 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
409 if(aodHandler->GetExtensions()) {
410
411 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
412 AliAODEvent *aodFromExt = ext->GetAOD();
413
414
415 if(fDecChannel==0){
416 absPdgMom=411;
417 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
418 }
419 if(fDecChannel==1){
420 absPdgMom=421;
421 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
422 }
423 if(fDecChannel==2){
424 absPdgMom=413;
425 arrayProng=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
426 }
427 }
428 } else {
429 if(fDecChannel==0){
430 absPdgMom=411;
431 arrayProng=(TClonesArray*)aod->GetList()->FindObject("Charm3Prong");
432 }
433 if(fDecChannel==1){
434 absPdgMom=421;
435 arrayProng=(TClonesArray*)aod->GetList()->FindObject("D0toKpi");
436 }
437 if(fDecChannel==2){
438 absPdgMom=413;
439 arrayProng=(TClonesArray*)aod->GetList()->FindObject("Dstar");
440 }
441 }
442
443 if(!arrayProng) {
444 AliError("AliAnalysisTaskSEHFv2::UserExec:Branch not found!\n");
445 return;
446 }
447
448 // fix for temporary bug in ESDfilter
449 // the AODs with null vertex pointer didn't pass the PhysSel
450 if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
451
452 TClonesArray *arrayMC=0;
453 AliAODMCHeader *mcHeader=0;
454
455 // load MC particles
456 if(fReadMC){
457
458 arrayMC = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
459 if(!arrayMC) {
460 AliWarning("AliAnalysisTaskSEHFv2::UserExec:MC particles branch not found!\n");
461 // return;
462 }
463
464 // load MC header
465 mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
466 if(!mcHeader) {
467 AliError("AliAnalysisTaskSEHFv2::UserExec:MC header branch not found!\n");
468 return;
469 }
470 }
471
472 fhEventsInfo->Fill(0); // count event
473
474 AliAODRecoDecayHF *d=0;
475
476 Int_t nCand = arrayProng->GetEntriesFast();
477 if(fDebug>2) printf("Number of D2H: %d\n",nCand);
478
479 // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
480 TString trigclass=aod->GetFiredTriggerClasses();
481 if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fhEventsInfo->Fill(5);
482
483 if(fRDCuts->IsEventSelectedInCentrality(aod)>0) return;
484 else fhEventsInfo->Fill(6);
485 if(fRDCuts->IsEventSelected(aod)) fhEventsInfo->Fill(1);
486 else{
487 if(fRDCuts->GetWhyRejection()==1) // rejected for pileup
488 fhEventsInfo->Fill(4);
489 return;
490 }
491
492 AliEventplane *pl=0x0;
493 TVector2* q=0x0;
494 Double_t rpangleevent=0;
495 Double_t eventplane=0;
496 //determine centrality bin
497 Float_t centr=fRDCuts->GetCentrality(aod);
498 Int_t icentr=0;
499 for(Int_t ic=5;ic<=80;ic=ic+5){
500 if(ic>centr){
501 icentr=ic;
502 break;
503 }
504 }
505 TString centrbinname=Form("centr%d_%d",icentr-5,icentr);
506
507 if(fUseV0EP){
508 rpangleevent=GetEventPlaneFromV0(aod);
509 eventplane=rpangleevent;
510 }else{
511 // event plane and resolution
512 //--------------------------------------------------------------------------
513 // extracting Q vectors and calculating v2 + resolution
514 pl = aod->GetHeader()->GetEventplaneP();
515 if(!pl){
516 AliError("AliAnalysisTaskSEHFv2::UserExec:no eventplane! v2 analysis without eventplane not possible!\n");
517 return;
518 }
519 q = pl->GetQVector();
520 rpangleevent = pl->GetEventplane("Q"); // reaction plane angle without autocorrelations removal
521 Double_t deltaPsi = pl->GetQsubRes();
522 if(TMath::Abs(deltaPsi)>TMath::Pi()/2.){
523 if(deltaPsi>0.) deltaPsi-=TMath::Pi();
524 else deltaPsi +=TMath::Pi();
525 } // difference of subevents reaction plane angle cannot be bigger than phi/2
526 Double_t planereso = TMath::Cos(2.*deltaPsi); // reaction plane resolution
527 //--------------------------------------------------------------------------
528 ((TH1F*)fOutput->FindObject(Form("hEvPlaneReso%s",centrbinname.Data())))->Fill(planereso);
529 }
530 ((TH1F*)fOutput->FindObject(Form("hEvPlane%s",centrbinname.Data())))->Fill(rpangleevent);
531
532 for (Int_t iCand = 0; iCand < nCand; iCand++) {
533
534 d=(AliAODRecoDecayHF*)arrayProng->UncheckedAt(iCand);
535
536 Bool_t isFidAcc = fRDCuts->IsInFiducialAcceptance(d->Pt(),d->Y(absPdgMom));
537 Int_t isSelected= fRDCuts->IsSelected(d,AliRDHFCuts::kCandidate,aod);
538 Bool_t isSelBit=kTRUE;
539 if(fDecChannel==0) isSelBit=d->HasSelectionBit(AliRDHFCuts::kDplusCuts);
540 if(fDecChannel==1) isSelBit=d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts);
541 if(fDecChannel==2) isSelBit=d->HasSelectionBit(AliRDHFCuts::kD0fromDstarCuts);
542 if(isSelected && isFidAcc && isSelBit) {
543 fhEventsInfo->Fill(2); // candidate selected
544 if(fDebug>3) printf("+++++++Is Selected\n");
545
546 Float_t* invMass=0x0;
547 Int_t nmasses;
548 CalculateInvMasses(d,invMass,nmasses);
549
550 Int_t ptbin=fRDCuts->PtBin(d->Pt());
551 if(ptbin==-1) {
552 fhEventsInfo->Fill(3);
553 continue;
554 }
555
556 if(!fUseV0EP) {
557 eventplane = GetEventPlaneForCandidate(d,q,pl); // remove autocorrelations
558 ((TH1F*)fOutput->FindObject(Form("hEvPlaneCand%s",centrbinname.Data())))->Fill(rpangleevent-eventplane);
559 //((TH1F*)fOutput->FindObject(Form("hEvPlaneCheck%s",centrbinname.Data())))->Fill((rpangleevent-eventplane)/100.*rpangleevent);
560 }
561
562 Float_t phi=d->Phi();
563 ((TH1F*)fOutput->FindObject(Form("hPhi_pt%d%s",ptbin,centrbinname.Data())))->Fill(phi);
564
565 Float_t deltaphi=GetPhi0Pi(phi-eventplane);
566
567 //fill the histograms with the appropriate method
568 if(fDecChannel==0)FillDplus(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr);
569 if(fDecChannel==1)FillD02p(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr);
570 if(fDecChannel==2)FillDstar(d,arrayMC,ptbin,deltaphi,invMass,isSelected,icentr);
571 }// end if selected
572 }
573
574 PostData(1,fhEventsInfo);
575 PostData(2,fOutput);
576 //PostData(4,flowEvent);
577 return;
578}
579
580//***************************************************************************
581
582// Methods used in the UserExec
583
584void AliAnalysisTaskSEHFv2::CalculateInvMasses(AliAODRecoDecayHF* d,Float_t*& masses,Int_t& nmasses){
585 //Calculates all the possible invariant masses for each candidate
586 //NB: the implementation for each candidate is responsibility of the corresponding developer
587
588 if(fDecChannel==0){
589 //D+ -- Giacomo
590 nmasses=1;
591 masses=new Float_t[nmasses];
592 Int_t pdgdaughters[3] = {211,321,211};
593 masses[0]=d->InvMass(3,(UInt_t*)pdgdaughters);
594 }
595 if(fDecChannel==1){
596 //D0 (Kpi) -- Chiara
597 const Int_t ndght=2;
598 nmasses=2;
599 masses=new Float_t[nmasses];
600 Int_t pdgdaughtersD0[ndght]={211,321};//pi,K
601 masses[0]=d->InvMass(ndght,(UInt_t*)pdgdaughtersD0); //D0
602 Int_t pdgdaughtersD0bar[ndght]={321,211};//K,pi
603 masses[1]=d->InvMass(ndght,(UInt_t*)pdgdaughtersD0bar); //D0bar
604 }
605 if(fDecChannel==2){
606 //D* -- Robert,Yifei, Alessandro
607 nmasses=1;
608 masses=new Float_t[nmasses];
609 masses[0]=((AliAODRecoCascadeHF*)d)->DeltaInvMass();
610 }
611}
612
613//******************************************************************************
614
615//Methods to fill the istograms with MC information, one for each candidate
616//NB: the implementation for each candidate is responsibility of the corresponding developer
617
618void AliAnalysisTaskSEHFv2::FillDplus(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi,Float_t* masses,Int_t isSel,Int_t icentr){
619 //D+ channel
620 if(!isSel){
621 if(fDebug>3)AliWarning("Candidate not selected\n");
622 return;
623 }
624 if(!masses){
625 if(fDebug>3)AliWarning("Masses not calculated\n");
626 return;
627 }
628
629 Int_t phibin=GetPhiBin(deltaphi);
630
631 ((TH1F*)fOutput->FindObject(Form("hMass_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
632 ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
633 ((TH2F*)fOutput->FindObject(Form("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[0]);
634 Int_t pdgdaughters[3] = {211,321,211};
635
636 if(fReadMC){
637 Int_t lab=-1;
638 lab = d->MatchToMC(411,arrayMC,3,pdgdaughters);
639 if(lab>=0){ //signal
640 ((TH1F*)fOutput->FindObject(Form("hSgn_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
641 ((TH2F*)fOutput->FindObject(Form("hMc2phiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
642 ((TH2F*)fOutput->FindObject(Form("hMphiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
643 } else{ //background
644 ((TH1F*)fOutput->FindObject(Form("hBkg_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
645 ((TH2F*)fOutput->FindObject(Form("hMc2phiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
646 ((TH2F*)fOutput->FindObject(Form("hMphiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
647 }
648 }
649}
650
651
652void AliAnalysisTaskSEHFv2::FillD02p(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi,Float_t* masses,Int_t isSel,Int_t icentr){
653
654 //D0->Kpi channel
655
656 //mass histograms
657 if(!masses){
658 if(fDebug>3)AliWarning("Masses not calculated\n");
659 return;
660 }
661 Int_t phibin=GetPhiBin(deltaphi);
662 if(isSel==1 || isSel==3) {
663 ((TH1F*)fOutput->FindObject(Form("hMass_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
664 ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
665 ((TH2F*)fOutput->FindObject(Form("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[0]);
666 }
667 if(isSel>=2) {
668 ((TH1F*)fOutput->FindObject(Form("hMass_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[1]);
669 ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
670 ((TH2F*)fOutput->FindObject(Form("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[1]);
671 }
672
673
674 //MC histograms
675 if(fReadMC){
676
677 Int_t matchtoMC=-1;
678
679 //D0
680 Int_t pdgdaughters[2];
681 pdgdaughters[0]=211;//pi
682 pdgdaughters[1]=321;//K
683 Int_t nprongs=2;
684 Int_t absPdgMom=421;
685
686 matchtoMC = d->MatchToMC(absPdgMom,arrayMC,nprongs,pdgdaughters);
687
688 Int_t prongPdgPlus=421,prongPdgMinus=(-1)*421;
689 if((isSel==1 || isSel==3)){ //D0
690 if(matchtoMC>=0){
691 AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
692 Int_t pdgMC = dMC->GetPdgCode();
693
694 if(pdgMC==prongPdgPlus) {
695 ((TH1F*)fOutput->FindObject(Form("hSgn_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
696 ((TH2F*)fOutput->FindObject(Form("hMc2phiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
697 ((TH2F*)fOutput->FindObject(Form("hMphiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
698 }
699 else {
700 ((TH1F*)fOutput->FindObject(Form("hRfl_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
701 ((TH2F*)fOutput->FindObject(Form("hMc2phiR_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
702 ((TH2F*)fOutput->FindObject(Form("hMphiRcentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[0]);
703 }
704 } else {
705 ((TH1F*)fOutput->FindObject(Form("hBkg_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
706 ((TH2F*)fOutput->FindObject(Form("hMc2phiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
707 ((TH2F*)fOutput->FindObject(Form("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[0]);
708 }
709 }
710 if(isSel>=2){ //D0bar
711 if(matchtoMC>=0){
712 AliAODMCParticle *dMC = (AliAODMCParticle*)arrayMC->At(matchtoMC);
713 Int_t pdgMC = dMC->GetPdgCode();
714
715 if(pdgMC==prongPdgMinus) {
716 ((TH1F*)fOutput->FindObject(Form("hSgn_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[1]);
717 ((TH2F*)fOutput->FindObject(Form("hMc2phiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
718 ((TH2F*)fOutput->FindObject(Form("hMphiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[1]);
719 }
720 else {
721 ((TH1F*)fOutput->FindObject(Form("hRfl_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[1]);
722 ((TH2F*)fOutput->FindObject(Form("hMc2phiR_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
723 ((TH2F*)fOutput->FindObject(Form("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[1]);
724 }
725 } else {
726 ((TH1F*)fOutput->FindObject(Form("hBkg_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[1]);
727 ((TH2F*)fOutput->FindObject(Form("hMc2phiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[1]);
728 ((TH2F*)fOutput->FindObject(Form("hMphiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[1]);
729 }
730 }
731 }
732}
733
734void AliAnalysisTaskSEHFv2::FillDstar(AliAODRecoDecayHF* d,TClonesArray *arrayMC,Int_t ptbin,Float_t deltaphi, Float_t* masses,Int_t isSel,Int_t icentr){
735 //D* channel
736 if(!isSel){
737 if(fDebug>3)AliWarning("Candidate not selected\n");
738 return;
739 }
740 if(!masses){
741 if(fDebug>3)AliWarning("Masses not calculated\n");
742 return;
743 }
744 Int_t phibin=GetPhiBin(deltaphi);
745
746 ((TH1F*)fOutput->FindObject(Form("hMass_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
747 ((TH2F*)fOutput->FindObject(Form("hMc2phi_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
748 ((TH2F*)fOutput->FindObject(Form("hMphicentr%d_%d",icentr-5,icentr)))->Fill(deltaphi,masses[0]);
749 Int_t pdgDgDStartoD0pi[2]={421,211};
750 Int_t pdgDgD0toKpi[2]={321,211};
751
752 if(fReadMC){
753 Int_t lab=-1;
754 lab = ((AliAODRecoCascadeHF*)d)->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,arrayMC);
755 if(lab>=0){ //signal
756 ((TH1F*)fOutput->FindObject(Form("hSgn_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
757 ((TH2F*)fOutput->FindObject(Form("hMc2phiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
758 ((TH2F*)fOutput->FindObject(Form("hMphiS_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
759 } else{ //background
760 ((TH1F*)fOutput->FindObject(Form("hBkg_pt%dphi%dcentr%d_%d",ptbin,phibin,icentr-5,icentr)))->Fill(masses[0]);
761 ((TH2F*)fOutput->FindObject(Form("hMc2phiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(TMath::Cos(2*deltaphi),masses[0]);
762 ((TH2F*)fOutput->FindObject(Form("hMphiB_pt%dcentr%d_%d",ptbin,icentr-5,icentr)))->Fill(deltaphi,masses[0]);
763 }
764 }
765
766
767}
768
769
770//________________________________________________________________________
771Int_t AliAnalysisTaskSEHFv2::GetPhiBin(Float_t deltaphi){
772
773 //give the bin corresponding to the value of deltaphi according to the binning requested in the constructor
774 Int_t phibin=0;
775 for(Int_t i=0;i<fNPhiBinLims-1;i++) {
776 if(deltaphi>=fPhiBins[i] && deltaphi<fPhiBins[i+1]) {
777 phibin=i;
778 break;
779 }
780 }
781 return phibin;
782}
783//________________________________________________________________________
784// Float_t AliAnalysisTaskSEHFv2::GetPhi02Pi(Float_t phi){
785// Float_t result=phi;
786// while(result<0){
787// result=result+2*TMath::Pi();
788// }
789// while(result>TMath::Pi()*2){
790// result=result-2*TMath::Pi();
791// }
792// return result;
793// }
794//________________________________________________________________________
795Float_t AliAnalysisTaskSEHFv2::GetPhi0Pi(Float_t phi){
796 Float_t result=phi;
797 while(result<0){
798 result=result+TMath::Pi();
799 }
800 while(result>TMath::Pi()){
801 result=result-TMath::Pi();
802 }
803 return result;
804}
805
806//________________________________________________________________________
807Float_t AliAnalysisTaskSEHFv2::GetEventPlaneForCandidate(AliAODRecoDecayHF* d, TVector2* q,AliEventplane *pl){
808 // remove autocorrelations
809
810 TArrayF* qx = pl->GetQContributionXArray();
811 TArrayF* qy = pl->GetQContributionYArray();
812 TVector2 qcopy = *q;
813
814 if(fDecChannel==2){
815 //D* -- Yifei, Alessandro,Robert
816 AliAODRecoDecayHF2Prong* theD0particle = ((AliAODRecoCascadeHF*)d)->Get2Prong();
817 AliAODTrack *track0 = (AliAODTrack*)theD0particle->GetDaughter(0);
818 AliAODTrack *track1 = (AliAODTrack*)theD0particle->GetDaughter(1);
819 AliAODTrack *track2 = ((AliAODRecoCascadeHF*)d)->GetBachelor();
820 // reduce global q vector
821
822 TVector2 q0;
823 if((track0->GetID()) < qx->fN){
824 q0.Set(qx->At(track0->GetID()),qy->At(track0->GetID()));}
825
826 TVector2 q1;
827 if((track1->GetID()) < qx->fN){
828 q1.Set(qx->At(track1->GetID()),qy->At(track1->GetID()));}
829
830 TVector2 q2;
831 if((track2->GetID()) < qx->fN){
832 q2.Set(qx->At(track2->GetID()),qy->At(track2->GetID()));}
833
834 qcopy = qcopy -(q0+q1+q2);
835
836 }
837
838 // reduce Q vector for D+ and D0
839
840 if(fDecChannel==1){
841 //D0 -- Chiara
842 AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
843 AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1);
844
845 TVector2 q0;
846 if((track0->GetID()) < qx->fN){
847 q0.Set(qx->At(track0->GetID()),qy->At(track0->GetID()));}
848
849 TVector2 q1;
850 if((track1->GetID()) < qx->fN){
851 q1.Set(qx->At(track1->GetID()),qy->At(track1->GetID()));}
852
853 qcopy = qcopy -(q0+q1);
854 }
855
856 if(fDecChannel==0){
857 //D+ -- Giacomo
858 AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
859 AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1);
860 AliAODTrack *track2 = (AliAODTrack*)d->GetDaughter(2);
861
862 TVector2 q0;
863 if((track0->GetID()) < qx->fN){
864 q0.Set(qx->At(track0->GetID()),qy->At(track0->GetID()));}
865
866 TVector2 q1;
867 if((track1->GetID()) < qx->fN){
868 q1.Set(qx->At(track1->GetID()),qy->At(track1->GetID()));}
869
870 TVector2 q2;
871 if((track2->GetID()) < qx->fN){
872 q2.Set(qx->At(track2->GetID()),qy->At(track2->GetID()));}
873
874 qcopy = qcopy -(q0+q1+q2);
875
876 }
877
878 return qcopy.Phi()/2.;
879
880}
881//________________________________________________________________________
882Float_t AliAnalysisTaskSEHFv2::GetEventPlaneFromV0(AliAODEvent *aodEvent){
883
884 Int_t centr=fRDCuts->GetCentrality(aodEvent);
885 centr=centr-centr%10;
886 //temporary fix
887 if(centr<20)centr=20;
888 if(centr>70)centr=70;
889 //end temporary fix
890 Int_t binx=0;
891 Int_t iParHist=(centr-20)/10;
892
893 TString name;name.Form("parhist%d_%d",centr,centr+10);
894
895 if(fDebug>15)printf("EPfromV0 centr %d, iparhist %d (%p-%p)\n",centr,iParHist,fParHist->FindObject(name.Data()),fParHist->At(iParHist));
896
897 Int_t runnumber=aodEvent->GetRunNumber();
898 if(fParHist->At(iParHist)){
899 for(Int_t i=1;i<=((TH2D*)fParHist->At(iParHist))->GetNbinsX()&&binx<=0;i++){
900 Int_t run=atoi(((TH2D*)fParHist->At(iParHist))->GetXaxis()->GetBinLabel(i));
901 if(run>=runnumber)binx=i;
902 }
903 }else{
904 fhEventsInfo->Fill(7);
905 }
906
907 AliFlowTrackCuts* cutsRP = AliFlowTrackCuts::GetStandardVZEROOnlyTrackCuts();
908 cutsRP->SetEvent(aodEvent, MCEvent());//, 0x0);
909 cutsRP->SetName( Form("rp_cuts") );
910 AliFlowTrackCuts* dummy = new AliFlowTrackCuts("null_cuts");
911 dummy->SetParamType(AliFlowTrackCuts::kGlobal);
912 dummy->SetPtRange(+1,-1); // select nothing QUICK
913 dummy->SetEtaRange(+1,-1); // select nothing VZERO
914 dummy->SetEvent(aodEvent,MCEvent());
915
916 //////////////// construct the flow event container ////////////
917 AliFlowEvent flowEvent(cutsRP,dummy);
918 flowEvent.SetReferenceMultiplicity( 64 );
919 for(Int_t i=0;i<64&&binx>0;i++){
920 AliFlowTrack *flowTrack=flowEvent.GetTrack(i);
921 Double_t inte=((TH2D*)fParHist->At(iParHist))->Integral(binx,binx,i+1,i+1);
922 if(inte>0)flowTrack->SetWeight(flowTrack->Weight()/inte);
923 }
924 if(fDebug>15)printf("EPfromV0 flow tracks weights done\n");
925
926 AliFlowVector qvec=flowEvent.GetQ(fV0EPorder);
927 Double_t angleEP=(1./(Double_t)fV0EPorder)*qvec.Phi();
928 if(fDebug>15)printf("EPfromV0 phi %f\n",angleEP);
929 return angleEP;
930}
931//________________________________________________________________________
932void AliAnalysisTaskSEHFv2::Terminate(Option_t */*option*/)
933{
934
935 // Terminate analysis
936 //
937 if(fDebug > 1) printf("AnalysisTaskSEHFv2: Terminate() \n");
938 /*
939 fhEventsInfo = dynamic_cast<TH1F*> (GetOutputData(1));
940 if(!fhEventsInfo){
941 printf("Error: hEventsInfo not available\n");
942 return;
943 }
944
945 fOutput = dynamic_cast<TList*> (GetOutputData(2));
946 if (!fOutput) {
947 printf("ERROR: fOutput not available\n");
948 return;
949 }
950 switch(fDecChannel){
951 case 0:
952 fRDCuts = dynamic_cast<AliRDHFCutsDplustoKpipi*> (GetOutputData(3));
953 break;
954 case 1:
955 fRDCuts = dynamic_cast<AliRDHFCutsD0toKpi*> (GetOutputData(3));
956 break;
957 case 2:
958 fRDCuts = dynamic_cast<AliRDHFCutsDStartoKpipi*> (GetOutputData(3));
959 break;
960 default:
961 break;
962 }
963 if (!fRDCuts) {
964 printf("ERROR: fRDCuts not available\n");
965 return;
966 }
967
968 // TCanvas* cvex=new TCanvas("example","Example Mass");
969 // cvex->cd();
970 // ((TH1F*)fOutput->FindObject("hMass_pt0phi0"))->Draw();
971 // TCanvas* cv2d=new TCanvas("2d","mass-cos2phi");
972 // cv2d->cd();
973 // ((TH2F*)fOutput->FindObject("hMc2phi"))->Draw("colz");
974 // TCanvas *cstat=new TCanvas("cstat","Stat");
975 // cstat->SetGridy();
976 // cstat->cd();
977 // fhEventsInfo->Draw("htext0");
978
979 // TMultiGraph *multig = new TMultiGraph();
980 if(fDecChannel==2)return;
981 TGraphErrors *g[fNPtBins];
982 TH1F *h = new TH1F("h","h",100,0.,1.);
983 TString hname;
984 TString gname;
985 for(Int_t ipt = 0;ipt<fNPtBins;ipt++){
986 g[ipt] = new TGraphErrors(fNPhiBinLims);
987 gname.Form("hMass_pt%d",ipt);
988 g[ipt]->SetTitle(gname.Data());
989 for(Int_t iphi = 0;iphi<fNPhiBinLims;iphi++){
990 hname.Form("hMass_pt%dphi%d",ipt,iphi);
991 h = (TH1F*)fOutput->FindObject("hMass_pt0phi0");
992 AliHFMassFitter fitter(h,fLowmasslimit,fUpmasslimit,2,0);
993 Int_t pdg=0;
994 switch(fDecChannel){
995 case 0:
996 pdg=411;
997 break;
998 case 1:
999 pdg=421;
1000 break;
1001 case 2:
1002 pdg=413;
1003 break;
1004 default:
1005 break;
1006 }
1007 fitter.SetInitialGaussianMean(TDatabasePDG::Instance()->GetParticle(pdg)->Mass());
1008 fitter.SetInitialGaussianSigma(0.012);
1009 fitter.InitNtuParam(Form("ntuPtbin%d",ipt));
1010 Double_t signal=0, errSignal=0;
1011 if(fitter.MassFitter(kFALSE)){
1012 fitter.Signal(3,signal,errSignal);
1013 }
1014 g[ipt]->SetPoint(iphi,fPhiBins[iphi],signal);
1015 g[ipt]->SetPointError(iphi,fPhiBins[iphi],errSignal);
1016 }//end loop over phi
1017 // multig->Add(g[ipt],"ap");
1018 }//end loop on pt
1019 TCanvas *cdndphi = new TCanvas("dN/d#phi","dN/d#phi");
1020 cdndphi->Divide(1,fNPtBins);
1021 for(Int_t ipt = 0;ipt<fNPtBins;ipt++){
1022 cdndphi->cd(ipt+1);
1023 g[ipt]->Draw("AP");
1024 }
1025 */
1026 return;
1027}
1028
1029//-------------------------------------------
1030/*
1031 Float_t GetEventPlaneFromVZERO(){
1032 AliAODHFUtil *info = (AliAODHFUtil*) aodEvent->GetList()->FindObject("fHFUtilInfo");
1033 if (!info) return -999.;
1034 //============= FIX ONLY FOR AOD033
1035 Double_t *par0;
1036 Double_t par0_137161[64] = { 6.71e-02 , 6.86e-02 , 7.06e-02 , 6.32e-02 ,
1037 5.91e-02 , 6.07e-02 , 5.78e-02 , 5.73e-02 , 5.91e-02 , 6.22e-02 ,
1038 5.90e-02 , 6.11e-02 , 5.55e-02 , 5.29e-02 , 5.19e-02 , 5.56e-02 ,
1039 6.25e-02 , 7.03e-02 , 5.64e-02 , 5.81e-02 , 4.57e-02 , 5.30e-02 ,
1040 5.13e-02 , 6.43e-02 , 6.27e-02 , 6.48e-02 , 6.07e-02 , 1.01e-01 ,
1041 6.68e-02 , 7.16e-02 , 6.36e-02 , 5.95e-02 , 2.52e-02 , 2.82e-02 ,
1042 2.56e-02 , 2.86e-02 , 2.82e-02 , 2.10e-02 , 2.13e-02 , 2.32e-02 ,
1043 2.75e-02 , 4.34e-02 , 3.78e-02 , 4.52e-02 , 4.11e-02 , 3.89e-02 ,
1044 4.10e-02 , 3.73e-02 , 4.51e-02 , 5.07e-02 , 5.42e-02 , 4.74e-02 ,
1045 4.33e-02 , 4.44e-02 , 4.64e-02 , 3.01e-02 , 6.38e-02 , 5.26e-02 ,
1046 4.99e-02 , 5.26e-02 , 5.47e-02 , 3.84e-02 , 5.00e-02 , 5.20e-02 };
1047 Double_t par0_137366[64] = { 7.12e-02 , 7.34e-02 , 7.39e-02 , 6.54e-02 , 6.11e-02 , 6.31e-02 , 6.15e-02 ,
1048 6.00e-02 , 6.10e-02 , 6.49e-02 , 6.17e-02 , 6.33e-02 , 6.00e-02 , 5.48e-02 ,
1049 5.44e-02 , 5.81e-02 , 6.49e-02 , 7.07e-02 , 5.91e-02 , 6.18e-02 , 4.82e-02 ,
1050 5.67e-02 , 5.36e-02 , 6.60e-02 , 6.37e-02 , 6.78e-02 , 6.31e-02 , 1.04e-01 ,
1051 6.91e-02 , 7.32e-02 , 6.61e-02 , 6.16e-02 , 2.64e-02 , 2.81e-02 , 2.64e-02 ,
1052 2.85e-02 , 2.87e-02 , 2.18e-02 , 2.19e-02 , 2.43e-02 , 2.81e-02 , 4.37e-02 ,
1053 3.90e-02 , 4.66e-02 , 4.24e-02 , 4.09e-02 , 4.21e-02 , 3.88e-02 , 4.83e-02 ,
1054 5.23e-02 , 5.44e-02 , 4.85e-02 , 4.42e-02 , 4.58e-02 , 4.74e-02 , 3.14e-02 ,
1055 6.31e-02 , 5.30e-02 , 5.01e-02 , 5.33e-02 , 5.70e-02 , 3.95e-02 , 4.98e-02 , 5.31e-02 };
1056 if (aodEvent->GetRunNumber() <= 137165)
1057 par0=par0_137161;
1058 else
1059 par0=par0_137366;
1060 Float_t vChCorr[64];
1061 for(int i=0; i!=64; ++i)
1062 vChCorr[i] = (info->GetVZEROChannel(i))/par0[i]/64.;
1063 //============= END OF FIX AOD033
1064 Float_t multR[8];
1065 for(int i=0; i!=8; ++i) multR[i] = 0;
1066 for(int i=0; i!=64; ++i)
1067 multR[i/8] += vChCorr[i];
1068 for(int i=0; i!=8; ++i)
1069 if(multR[i]) {
1070 double Qx=0, Qy=0;
1071 for(int j=0; j!=8; ++j) {
1072 double phi = TMath::PiOver4()*(0.5+j);
1073 Qx += TMath::Cos(2*phi)*vChCorr[i*8+j]/multR[i];
1074 Qy += TMath::Sin(2*phi)*vChCorr[i*8+j]/multR[i];
1075 }
1076 }
1077 return 0.5*TMath::ATan2(Qy,Qx)+TMath::PiOver2();
1078 }
1079
1080*/