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