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