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