]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliCFTaskVertexingHF.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Including: 1) external setter of MC option and proper usage on the task, 2) 3D histog...
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliCFTaskVertexingHF.cxx
CommitLineData
379592af 1/**************************************************************************
2 * Copyright(c) 1998-2009, 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// Class for HF corrections as a function of many variables
18// 6 Steps introduced: MC, MC Acc, Reco, Reco Acc, Reco Acc + ITS Cl,
19// Reco Acc + ITS Cl + PPR cuts
20// 12 variables used: pt, y, cosThetaStar, ptPi, ptK, ct,
21// dca, d0Pi, d0K, d0Pixd0K, cosPointingAngle, phi
22//
23//-----------------------------------------------------------------------
24// Author : C. Zampolli, CERN
25// D. Caffarri, Univ & INFN Padova caffarri@pd.infn.it
26//-----------------------------------------------------------------------
27//-----------------------------------------------------------------------
28// Base class for HF Unfolding (pt and eta)
29// correlation matrix filled at Acceptance and PPR level
30// Author: A.Grelli , Utrecht - agrelli@uu.nl
31//-----------------------------------------------------------------------
32#include <TCanvas.h>
33#include <TParticle.h>
34#include <TDatabasePDG.h>
35#include <TH1I.h>
36#include <TStyle.h>
37#include <TFile.h>
5b37c6e5 38#include <TF1.h>
379592af 39
40#include "AliCFTaskVertexingHF.h"
41#include "AliStack.h"
42#include "AliMCEvent.h"
43#include "AliCFManager.h"
44#include "AliCFContainer.h"
45#include "AliLog.h"
5cd139bc 46#include "AliInputEventHandler.h"
379592af 47#include "AliAnalysisManager.h"
48#include "AliAODHandler.h"
49#include "AliAODEvent.h"
50#include "AliAODRecoDecay.h"
51#include "AliAODRecoDecayHF.h"
52#include "AliAODRecoDecayHF2Prong.h"
f2703bd2 53#include "AliAODRecoDecayHF3Prong.h"
54#include "AliAODRecoDecayHF4Prong.h"
55#include "AliAODRecoCascadeHF.h"
379592af 56#include "AliAODMCParticle.h"
57#include "AliAODMCHeader.h"
58#include "AliESDtrack.h"
59#include "TChain.h"
60#include "THnSparse.h"
61#include "TH2D.h"
62#include "AliESDtrackCuts.h"
f2703bd2 63#include "AliRDHFCuts.h"
379592af 64#include "AliRDHFCutsD0toKpi.h"
f2703bd2 65#include "AliRDHFCutsDplustoKpipi.h"
66#include "AliRDHFCutsDStartoKpipi.h"
67#include "AliRDHFCutsDstoKKpi.h"
68#include "AliRDHFCutsLctopKpi.h"
69#include "AliRDHFCutsD0toKpipipi.h"
5cd139bc 70#include "AliRDHFCutsLctoV0.h"
379592af 71#include "AliCFVertexingHF2Prong.h"
4e600f58 72#include "AliCFVertexingHF3Prong.h"
367e9aa3 73#include "AliCFVertexingHFCascade.h"
5cd139bc 74#include "AliCFVertexingHFLctoV0bachelor.h"
379592af 75#include "AliCFVertexingHF.h"
afb24c40 76#include "AliVertexingHFUtils.h"
f2703bd2 77#include "AliAnalysisDataSlot.h"
78#include "AliAnalysisDataContainer.h"
5cd139bc 79#include "AliPIDResponse.h"
379592af 80
81//__________________________________________________________________________
82AliCFTaskVertexingHF::AliCFTaskVertexingHF() :
5cd139bc 83 AliAnalysisTaskSE(),
84 fCFManager(0x0),
85 fHistEventsProcessed(0x0),
86 fCorrelation(0x0),
87 fCountMC(0),
88 fCountAcc(0),
89 fCountVertex(0),
90 fCountRefit(0),
91 fCountReco(0),
92 fCountRecoAcc(0),
93 fCountRecoITSClusters(0),
94 fCountRecoPPR(0),
95 fCountRecoPID(0),
96 fEvents(0),
97 fDecayChannel(0),
98 fFillFromGenerated(kFALSE),
99 fOriginDselection(0),
100 fAcceptanceUnf(kTRUE),
101 fCuts(0),
102 fUseWeight(kFALSE),
103 fWeight(1.),
104 fUseFlatPtWeight(kFALSE),
105 fUseZWeight(kFALSE),
106 fUseNchWeight(kFALSE),
107 fNvar(0),
108 fPartName(""),
109 fDauNames(""),
110 fSign(2),
111 fCentralitySelection(kTRUE),
112 fFakeSelection(0),
113 fRejectIfNoQuark(kTRUE),
114 fUseMCVertex(kFALSE),
115 fDsOption(1),
116 fGenDsOption(3),
117 fConfiguration(kCheetah), // by default, setting the fast configuration
118 fFuncWeight(0x0),
119 fHistoMeasNch(0x0),
120 fHistoMCNch(0x0),
121 fResonantDecay(0),
122 fLctoV0bachelorOption(1),
ca7d039b 123 fGenLctoV0bachelorOption(0),
04b11915 124 fUseSelectionBit(kTRUE),
125 fPDGcode(0)
379592af 126{
5cd139bc 127 //
128 //Default ctor
129 //
379592af 130}
131//___________________________________________________________________________
5b37c6e5 132AliCFTaskVertexingHF::AliCFTaskVertexingHF(const Char_t* name, AliRDHFCuts* cuts, TF1* func) :
5cd139bc 133 AliAnalysisTaskSE(name),
134 fCFManager(0x0),
135 fHistEventsProcessed(0x0),
136 fCorrelation(0x0),
137 fCountMC(0),
138 fCountAcc(0),
139 fCountVertex(0),
140 fCountRefit(0),
141 fCountReco(0),
142 fCountRecoAcc(0),
143 fCountRecoITSClusters(0),
144 fCountRecoPPR(0),
145 fCountRecoPID(0),
146 fEvents(0),
147 fDecayChannel(0),
148 fFillFromGenerated(kFALSE),
149 fOriginDselection(0),
150 fAcceptanceUnf(kTRUE),
151 fCuts(cuts),
152 fUseWeight(kFALSE),
153 fWeight(1.),
154 fUseFlatPtWeight(kFALSE),
155 fUseZWeight(kFALSE),
156 fUseNchWeight(kFALSE),
157 fNvar(0),
158 fPartName(""),
159 fDauNames(""),
160 fSign(2),
161 fCentralitySelection(kTRUE),
162 fFakeSelection(0),
163 fRejectIfNoQuark(kTRUE),
164 fUseMCVertex(kFALSE),
165 fDsOption(1),
166 fGenDsOption(3),
167 fConfiguration(kCheetah), // by default, setting the fast configuration
168 fFuncWeight(func),
169 fHistoMeasNch(0x0),
170 fHistoMCNch(0x0),
171 fResonantDecay(0),
172 fLctoV0bachelorOption(1),
ca7d039b 173 fGenLctoV0bachelorOption(0),
04b11915 174 fUseSelectionBit(kTRUE),
175 fPDGcode(0)
379592af 176{
5cd139bc 177 //
178 // Constructor. Initialization of Inputs and Outputs
179 //
180 /*
181 DefineInput(0) and DefineOutput(0)
182 are taken care of by AliAnalysisTaskSE constructor
183 */
184 DefineOutput(1,TH1I::Class());
185 DefineOutput(2,AliCFContainer::Class());
186 DefineOutput(3,THnSparseD::Class());
187 DefineOutput(4,AliRDHFCuts::Class());
379592af 188
5cd139bc 189 fCuts->PrintAll();
379592af 190}
191
192//___________________________________________________________________________
193AliCFTaskVertexingHF& AliCFTaskVertexingHF::operator=(const AliCFTaskVertexingHF& c)
194{
5cd139bc 195 //
196 // Assignment operator
197 //
198 if (this!=&c) {
199 AliAnalysisTaskSE::operator=(c) ;
200 fCFManager = c.fCFManager;
201 fHistEventsProcessed = c.fHistEventsProcessed;
202 fCuts = c.fCuts;
203 fFuncWeight = c.fFuncWeight;
204 fHistoMeasNch = c.fHistoMeasNch;
205 fHistoMCNch = c.fHistoMCNch;
206 }
207 return *this;
379592af 208}
209
210//___________________________________________________________________________
211AliCFTaskVertexingHF::AliCFTaskVertexingHF(const AliCFTaskVertexingHF& c) :
5cd139bc 212 AliAnalysisTaskSE(c),
213 fCFManager(c.fCFManager),
214 fHistEventsProcessed(c.fHistEventsProcessed),
215 fCorrelation(c.fCorrelation),
216 fCountMC(c.fCountMC),
217 fCountAcc(c.fCountAcc),
218 fCountVertex(c.fCountVertex),
219 fCountRefit(c.fCountRefit),
220 fCountReco(c.fCountReco),
221 fCountRecoAcc(c.fCountRecoAcc),
222 fCountRecoITSClusters(c.fCountRecoITSClusters),
223 fCountRecoPPR(c.fCountRecoPPR),
224 fCountRecoPID(c.fCountRecoPID),
225 fEvents(c.fEvents),
226 fDecayChannel(c.fDecayChannel),
227 fFillFromGenerated(c.fFillFromGenerated),
228 fOriginDselection(c.fOriginDselection),
229 fAcceptanceUnf(c.fAcceptanceUnf),
230 fCuts(c.fCuts),
231 fUseWeight(c.fUseWeight),
232 fWeight(c.fWeight),
233 fUseFlatPtWeight(c.fUseFlatPtWeight),
234 fUseZWeight(c.fUseZWeight),
235 fUseNchWeight(c.fUseNchWeight),
236 fNvar(c.fNvar),
237 fPartName(c.fPartName),
238 fDauNames(c.fDauNames),
239 fSign(c.fSign),
240 fCentralitySelection(c.fCentralitySelection),
241 fFakeSelection(c.fFakeSelection),
242 fRejectIfNoQuark(c.fRejectIfNoQuark),
243 fUseMCVertex(c.fUseMCVertex),
244 fDsOption(c.fDsOption),
245 fGenDsOption(c.fGenDsOption),
246 fConfiguration(c.fConfiguration),
247 fFuncWeight(c.fFuncWeight),
248 fHistoMeasNch(c.fHistoMeasNch),
249 fHistoMCNch(c.fHistoMCNch),
250 fResonantDecay(c.fResonantDecay),
251 fLctoV0bachelorOption(c.fLctoV0bachelorOption),
ca7d039b 252 fGenLctoV0bachelorOption(c.fGenLctoV0bachelorOption),
04b11915 253 fUseSelectionBit(c.fUseSelectionBit),
254 fPDGcode(c.fPDGcode)
379592af 255{
5cd139bc 256 //
257 // Copy Constructor
258 //
379592af 259}
260
261//___________________________________________________________________________
f2703bd2 262AliCFTaskVertexingHF::~AliCFTaskVertexingHF()
263{
5cd139bc 264 //
265 //destructor
266 //
267 if (fCFManager) delete fCFManager ;
268 if (fHistEventsProcessed) delete fHistEventsProcessed ;
269 if (fCorrelation) delete fCorrelation ;
270 if (fCuts) delete fCuts;
271 if (fFuncWeight) delete fFuncWeight;
272 if (fHistoMeasNch) delete fHistoMeasNch;
273 if (fHistoMCNch) delete fHistoMCNch;
379592af 274}
275
f2703bd2 276//_________________________________________________________________________-
277void AliCFTaskVertexingHF::Init()
278{
5cd139bc 279 //
280 // Initialization
281 //
f2703bd2 282
5cd139bc 283 if (fDebug>1) printf("AliCFTaskVertexingHF::Init()");
284 if(fUseWeight && fUseZWeight) { AliFatal("Can not use at the same time pt and z-vtx weights, please choose"); return; }
285 if(fUseWeight && fUseNchWeight) { AliFatal("Can not use at the same time pt and Nch weights, please choose"); return; }
286 if(fUseNchWeight && !fHistoMCNch) { AliFatal("Need to pass the MC Nch distribution to use Nch weights"); return; }
287 if(fUseNchWeight) CreateMeasuredNchHisto();
8a983038 288
5cd139bc 289 AliRDHFCuts *copyfCuts = 0x0;
290 if (!fCuts){
291 AliFatal("No cuts defined - Exiting...");
292 return;
293 }
294
295 switch (fDecayChannel){
296 case 2:{
04b11915 297 fPDGcode = 421;
5cd139bc 298 copyfCuts = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
299 switch (fConfiguration) {
300 case kSnail: // slow configuration: all variables in
301 fNvar = 16;
302 break;
303 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
304 fNvar = 8;
305 break;
306 }
307 fPartName="D0";
308 fDauNames="K+pi";
309 break;
310 }
311 case 21:{
04b11915 312 fPDGcode = 413;
5cd139bc 313 copyfCuts = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
314 switch (fConfiguration) {
315 case kSnail: // slow configuration: all variables in
316 fNvar = 16;
317 break;
318 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
319 fNvar = 8;
320 break;
321 }
322 fPartName="Dstar";
323 fDauNames="K+pi+pi";
324 break;
325 }
326 case 22:{
04b11915 327 fPDGcode = 4122;
5cd139bc 328 copyfCuts = new AliRDHFCutsLctoV0(*(static_cast<AliRDHFCutsLctoV0*>(fCuts)));
329 switch (fConfiguration) {
330 case kSnail: // slow configuration: all variables in
331 fNvar = 16;
332 break;
333 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
334 fNvar = 8;
335 break;
336 }
337 fPartName="Lambdac";
338 fDauNames="V0+bachelor";
339 break;
340 }
341 case 31:{
04b11915 342 fPDGcode = 411;
5cd139bc 343 copyfCuts = new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
344 switch (fConfiguration) {
345 case kSnail: // slow configuration: all variables in
346 fNvar = 14;
347 break;
348 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
349 fNvar = 8;
350 break;
351 }
352 fPartName="Dplus";
353 fDauNames="K+pi+pi";
354 break;
355 }
356 case 32:{
04b11915 357 fPDGcode = 4122;
5cd139bc 358 copyfCuts = new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fCuts)));
359 switch (fConfiguration) {
360 case kSnail: // slow configuration: all variables in
361 fNvar = 18;
362 break;
363 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
364 fNvar = 8;
365 break;
366 }
367 fPartName="Lambdac";
368 fDauNames="p+K+pi";
369 break;
370 }
371 case 33:{
04b11915 372 fPDGcode = 431;
5cd139bc 373 copyfCuts = new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
374 switch (fConfiguration) {
375 case kSnail: // slow configuration: all variables in
376 fNvar = 14;
377 break;
378 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
379 fNvar = 8;
380 break;
381 }
382 fPartName="Ds";
383 fDauNames="K+K+pi";
384 break;
385 }
386 case 4:{
04b11915 387 fPDGcode = 421;
5cd139bc 388 copyfCuts = new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
389 switch (fConfiguration) {
390 case kSnail: // slow configuration: all variables in
391 fNvar = 16;
392 break;
393 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
394 fNvar = 8;
395 break;
396 }
397 fPartName="D0";
398 fDauNames="K+pi+pi+pi";
399 break;
400 }
401 default:
402 AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exiting...");
403 break;
404 }
f2703bd2 405
5cd139bc 406 const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
407 if (copyfCuts){
408 copyfCuts->SetName(nameoutput);
8a983038 409
5cd139bc 410 //Post the data
411 PostData(4, copyfCuts);
412 }
413 else{
414 AliFatal("Failing initializing AliRDHFCuts object - Exiting...");
415 }
f2703bd2 416
5cd139bc 417 return;
f2703bd2 418}
419
379592af 420//_________________________________________________
421void AliCFTaskVertexingHF::UserExec(Option_t *)
422{
5cd139bc 423 //
424 // Main loop function
425 //
379592af 426
5cd139bc 427 PostData(1,fHistEventsProcessed) ;
428 PostData(2,fCFManager->GetParticleContainer()) ;
429 PostData(3,fCorrelation) ;
b7af2639 430
5cd139bc 431 if (fFillFromGenerated){
432 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
433 }
379592af 434
5cd139bc 435 if (!fInputEvent) {
436 Error("UserExec","NO EVENT FOUND!");
437 return;
438 }
379592af 439
5cd139bc 440 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
379592af 441
5cd139bc 442 TClonesArray *arrayBranch=0;
379592af 443
5cd139bc 444 if(!aodEvent && AODEvent() && IsStandardAOD()) {
445 // In case there is an AOD handler writing a standard AOD, use the AOD
446 // event in memory rather than the input (ESD) event.
447 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
448 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
449 // have to taken from the AOD event hold by the AliAODExtension
450 AliAODHandler* aodHandler = (AliAODHandler*)
451 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
452 if(aodHandler->GetExtensions()) {
453 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
454 AliAODEvent *aodFromExt = ext->GetAOD();
f2703bd2 455
5cd139bc 456 switch (fDecayChannel){
457 case 2:{
458 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
459 break;
460 }
461 case 21:{
462 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
463 break;
464 }
465 case 22:{
466 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("CascadesHF");
467 break;
468 }
469 case 31:
470 case 32:
471 case 33:{
472 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
473 break;
474 }
475 case 4:{
476 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
477 break;
478 }
479 default:
480 break;
481 }
482 }
483 }
484 else {
485 switch (fDecayChannel){
486 case 2:{
487 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
488 break;
489 }
490 case 21:{
491 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
492 break;
493 }
494 case 22:{
495 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("CascadesHF");
496 break;
497 }
498 case 31:
499 case 32:
500 case 33:{
501 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
502 break;
503 }
504 case 4:{
505 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm4Prong");
506 break;
507 }
508 default:
509 break;
510 }
511 }
ec5288c3 512
5cd139bc 513 AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
514 if (!aodVtx) return;
f2703bd2 515
5cd139bc 516 if (!arrayBranch) {
517 AliError("Could not find array of HF vertices");
518 return;
519 }
379592af 520
5cd139bc 521 fEvents++;
6e2e4f50 522
5cd139bc 523 fCFManager->SetRecEventInfo(aodEvent);
524 fCFManager->SetMCEventInfo(aodEvent);
379592af 525
5cd139bc 526 //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro.
379592af 527
5cd139bc 528 //loop on the MC event
379592af 529
5cd139bc 530 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
531 if (!mcArray) {
532 AliError("Could not find Monte-Carlo in AOD");
533 return;
534 }
535 Int_t icountMC = 0;
536 Int_t icountAcc = 0;
537 Int_t icountReco = 0;
538 Int_t icountVertex = 0;
539 Int_t icountRefit = 0;
540 Int_t icountRecoAcc = 0;
541 Int_t icountRecoITSClusters = 0;
542 Int_t icountRecoPPR = 0;
543 Int_t icountRecoPID = 0;
544 Int_t cquarks = 0;
f2703bd2 545
5cd139bc 546 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
547 if (!mcHeader) {
548 AliError("Could not find MC Header in AOD");
549 return;
550 }
e11ae259 551
78471b33 552 fHistEventsProcessed->Fill(0.5);
553
5cd139bc 554 Double_t* containerInput = new Double_t[fNvar];
555 Double_t* containerInputMC = new Double_t[fNvar];
e11ae259 556
f2703bd2 557
5cd139bc 558 AliCFVertexingHF* cfVtxHF=0x0;
559 switch (fDecayChannel){
560 case 2:{
561 cfVtxHF = new AliCFVertexingHF2Prong(mcArray, fOriginDselection);
562 break;
563 }
564 case 21:{
565 cfVtxHF = new AliCFVertexingHFCascade(mcArray, fOriginDselection);
566 break;
567 }
568 case 22:{
569 cfVtxHF = new AliCFVertexingHFLctoV0bachelor(mcArray, fOriginDselection,fGenLctoV0bachelorOption); // Lc -> K0S+proton
570 break;
571 }
572 case 31:
573 // case 32:
574 case 33:{
575 cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel);
576 if(fDecayChannel==33){
577 ((AliCFVertexingHF3Prong*)cfVtxHF)->SetGeneratedDsOption(fGenDsOption);
578 }
579 break;
580 }
581 case 32:{
582 cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel,fResonantDecay);
583 }
584 case 4:{
585 //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection); // not there yet
586 break;
587 }
588 default:
589 break;
590 }
591 if (!cfVtxHF){
592 AliError("No AliCFVertexingHF initialized");
593 delete[] containerInput;
594 delete[] containerInputMC;
595 return;
596 }
f2703bd2 597
5cd139bc 598 Double_t zPrimVertex = aodVtx ->GetZ();
599 Double_t zMCVertex = mcHeader->GetVtxZ();
600 Int_t runnumber = aodEvent->GetRunNumber();
601 fWeight=1.;
602 if(fUseZWeight) fWeight *= GetZWeight(zMCVertex,runnumber);
603 if(fUseNchWeight){
604 Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(mcArray,-1.0,1.0);
605 fWeight *= GetNchWeight(nChargedMCPhysicalPrimary);
606 AliDebug(2,Form("Using Nch weights, Mult=%d Weight=%f\n",nChargedMCPhysicalPrimary,fWeight));
607 }
17bf1a34 608
5cd139bc 609 if (TMath::Abs(zMCVertex) > fCuts->GetMaxVtxZ()){
610 AliDebug(3,Form("z coordinate of MC vertex = %f, it was required to be within [-%f, +%f], skipping event", zMCVertex, fCuts->GetMaxVtxZ(), fCuts->GetMaxVtxZ()));
611 delete[] containerInput;
612 delete[] containerInputMC;
613 return;
614 }
751b8274 615
5cd139bc 616 AliESDtrackCuts** trackCuts = new AliESDtrackCuts*[cfVtxHF->GetNProngs()];
617 if (fDecayChannel == 21){
618 // for the D*, setting the third element of the array of the track cuts to those for the soft pion
619 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs()-1; iProng++){
620 trackCuts[iProng]=fCuts->GetTrackCuts();
621 }
622 trackCuts[2] = fCuts->GetTrackCutsSoftPi();
623 }
624 else if (fDecayChannel == 22) {
625 // for the Lc->V0+bachelor, setting the second and third elements of the array of the track cuts to those for the V0 daughters
626 trackCuts[0]=fCuts->GetTrackCuts();
627 trackCuts[1]=fCuts->GetTrackCutsV0daughters();
628 trackCuts[2]=fCuts->GetTrackCutsV0daughters();
629 }
630 else {
631 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs(); iProng++){
632 trackCuts[iProng]=fCuts->GetTrackCuts();
633 }
634 }
2bf2e62b 635
5cd139bc 636 //General settings: vertex, feed down and fill reco container with generated values.
637 cfVtxHF->SetRecoPrimVertex(zPrimVertex);
638 cfVtxHF->SetMCPrimaryVertex(zMCVertex);
639 cfVtxHF->SetFillFromGenerated(fFillFromGenerated);
640 cfVtxHF->SetNVar(fNvar);
641 cfVtxHF->SetFakeSelection(fFakeSelection);
642 cfVtxHF->SetRejectCandidateIfNotFromQuark(fRejectIfNoQuark);
643 cfVtxHF->SetConfiguration(fConfiguration);
644
645 // switch-off the trigger class selection (doesn't work for MC)
646 fCuts->SetTriggerClass("");
647
648 // MC vertex, to be used, in case, for pp
649 if (fUseMCVertex) fCuts->SetUseMCVertex();
650
651 if (fCentralitySelection){ // keep only the requested centrality
652 if(fCuts->IsEventSelectedInCentrality(aodEvent)!=0) {
653 delete[] containerInput;
654 delete[] containerInputMC;
655 delete [] trackCuts;
656 return;
657 }
658 } else { // keep all centralities
659 fCuts->SetMinCentrality(0.);
660 fCuts->SetMaxCentrality(100.);
661 }
8b1b055d 662
5cd139bc 663 Float_t centValue = fCuts->GetCentrality(aodEvent);
664 cfVtxHF->SetCentralityValue(centValue);
f2703bd2 665
5cd139bc 666 // number of tracklets - multiplicity estimator
bceabe60 667 Double_t multiplicity = (Double_t)(AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-1.,1.)); // casted to double because the CF is filled with doubles
5cd139bc 668 cfVtxHF->SetMultiplicity(multiplicity);
ec5288c3 669
5cd139bc 670 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
671 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
672 if (!mcPart){
673 AliError("Failed casting particle from MC array!, Skipping particle");
674 continue;
675 }
04b11915 676
677 //counting c quarks
678 cquarks += cfVtxHF->MCcquarkCounting(mcPart);
679
5cd139bc 680 // check the MC-level cuts, must be the desidered particle
681 if (!fCFManager->CheckParticleCuts(0, mcPart)) {
682 AliDebug(2,"Check the MC-level cuts - not desidered particle");
683 continue; // 0 stands for MC level
684 }
685 cfVtxHF->SetMCCandidateParam(iPart);
b7af2639 686
b7af2639 687
5cd139bc 688 if (!(cfVtxHF->SetLabelArray())){
689 AliDebug(2,Form("Impossible to set the label array (decaychannel = %d)",fDecayChannel));
690 continue;
691 }
2bf2e62b 692
5cd139bc 693 //check the candiate family at MC level
694 if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) {
695 AliDebug(2,Form("Check on the family wrong!!! (decaychannel = %d)",fDecayChannel));
696 continue;
697 }
698 else{
699 AliInfo(Form("Check on the family OK!!! (decaychannel = %d)",fDecayChannel));
700 }
f2703bd2 701
5cd139bc 702 //Fill the MC container
703 Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
704 AliDebug(2,Form("mcContainerFilled = %d)",mcContainerFilled));
705 if (mcContainerFilled) {
706 if (fUseWeight){
707 if (fFuncWeight){ // using user-defined function
708 AliDebug(2,"Using function");
709 fWeight = fFuncWeight->Eval(containerInputMC[0]);
710 }
711 else{ // using FONLL
712 AliDebug(2,"Using FONLL");
713 fWeight = GetWeight(containerInputMC[0]);
714 }
715 AliDebug(2,Form("pt = %f, weight = %f",containerInputMC[0], fWeight));
716 }
717 if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) continue;
718 //MC Limited Acceptance
719 if (TMath::Abs(containerInputMC[1]) < 0.5) {
720 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight);
721 AliDebug(3,"MC Lim Acc container filled\n");
722 }
f2703bd2 723
5cd139bc 724 //MC
725 fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight);
726 icountMC++;
727 AliDebug(3,"MC cointainer filled \n");
f2703bd2 728
5cd139bc 729 // MC in acceptance
730 // check the MC-Acceptance level cuts
731 // since standard CF functions are not applicable, using Kine Cuts on daughters
732 Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep();
733 if (mcAccepStep){
734 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight);
735 AliDebug(3,"MC acceptance cut passed\n");
736 icountAcc++;
f2703bd2 737
5cd139bc 738 //MC Vertex step
739 if (fCuts->IsEventSelected(aodEvent)){
740 // filling the container if the vertex is ok
741 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex, fWeight) ;
742 AliDebug(3,"Vertex cut passed and container filled\n");
743 icountVertex++;
f2703bd2 744
5cd139bc 745 //mc Refit requirement
746 Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, &trackCuts[0]);
747 if (mcRefitStep){
748 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight);
749 AliDebug(3,"MC Refit cut passed and container filled\n");
750 icountRefit++;
751 }
752 else{
753 AliDebug(3,"MC Refit cut not passed\n");
754 continue;
755 }
379592af 756 }
5cd139bc 757 else{
758 AliDebug (3, "MC vertex step not passed\n");
759 continue;
760 }
761 }
762 else{
763 AliDebug (3, "MC in acceptance step not passed\n");
764 continue;
765 }
766 }
767 else {
768 AliDebug (3, "MC container not filled\n");
769 }
770 }
f2703bd2 771
5cd139bc 772 if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks));
773 AliDebug(2,Form("Found %i MC particles that are %s!!",icountMC,fPartName.Data()));
774 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Acc cuts!!",icountAcc,fPartName.Data()));
775 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Vertex cuts!!",icountVertex,fPartName.Data()));
776 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Refit cuts!!",icountRefit,fPartName.Data()));
777
778 // Now go to rec level
779 fCountMC += icountMC;
780 fCountAcc += icountAcc;
781 fCountVertex+= icountVertex;
782 fCountRefit+= icountRefit;
783
784 AliDebug(2,Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel));
379592af 785
5cd139bc 786 for(Int_t iCandid = 0; iCandid<arrayBranch->GetEntriesFast();iCandid++){
787 AliAODRecoDecayHF* charmCandidate=0x0;
788 switch (fDecayChannel){
789 case 2:{
790 charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid);
791 break;
792 }
793 case 21:
794 case 22:{
795 charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid);
796 break;
797 }
798 case 31:
799 case 32:
800 case 33:{
801 charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid);
802 break;
803 }
804 case 4:{
805 charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid);
806 break;
807 }
808 default:
809 break;
810 }
f2703bd2 811
5cd139bc 812 Bool_t unsetvtx=kFALSE;
813 if(!charmCandidate->GetOwnPrimaryVtx()) {
814 charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
815 unsetvtx=kTRUE;
816 }
f2703bd2 817
5cd139bc 818 Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate);
819 if (!signAssociation){
820 charmCandidate = 0x0;
821 continue;
822 }
3ee5eb83 823
5cd139bc 824 Int_t isPartOrAntipart = cfVtxHF->CheckReflexion(fSign);
825 if (isPartOrAntipart == 0){
826 AliDebug(2, Form("The candidate pdg code doesn't match the requirement set in the task (fSign = %d)",fSign));
827 continue;
828 }
3ee5eb83 829
5cd139bc 830 AliDebug(3,Form("iCandid=%d - signAssociation=%d, isPartOrAntipart=%d",iCandid, signAssociation, isPartOrAntipart));
6694a2a8 831
5cd139bc 832 Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
833 if (recoContFilled){
6694a2a8 834
5cd139bc 835 // weight according to pt
836 if (fUseWeight){
837 if (fFuncWeight){ // using user-defined function
838 AliDebug(2, "Using function");
839 fWeight = fFuncWeight->Eval(containerInput[0]);
840 }
841 else{ // using FONLL
842 AliDebug(2, "Using FONLL");
843 fWeight = GetWeight(containerInput[0]);
844 }
845 AliDebug(2, Form("pt = %f, weight = %f",containerInput[0], fWeight));
846 }
6694a2a8 847
5cd139bc 848 if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue;
f2703bd2 849
5cd139bc 850 //Reco Step
851 Bool_t recoStep = cfVtxHF->RecoStep();
852 Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
f2703bd2 853
5cd139bc 854
ca7d039b 855 // Selection on the filtering bit
856 Bool_t isBitSelected = kTRUE;
857 if(fDecayChannel==2) {
858 if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)) isBitSelected = kFALSE;
859 }else if(fDecayChannel==31){
860 if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kDplusCuts)) isBitSelected = kFALSE;
861 }else if(fDecayChannel==33){
862 if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kDsCuts)) isBitSelected = kFALSE;
863 }
864 if(!isBitSelected) continue;
865
866
867
5cd139bc 868 if (recoStep && recoContFilled && vtxCheck){
869 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
870 icountReco++;
871 AliDebug(3,"Reco step passed and container filled\n");
f2703bd2 872
5cd139bc 873 //Reco in the acceptance -- take care of UNFOLDING!!!!
874 Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(&trackCuts[0]);
875 if (recoAcceptanceStep) {
876 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
877 icountRecoAcc++;
878 AliDebug(3,"Reco acceptance cut passed and container filled\n");
f2703bd2 879
5cd139bc 880 if(fAcceptanceUnf){
881 Double_t fill[4]; //fill response matrix
04b11915 882 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fPDGcode,fill);
5cd139bc 883 if (bUnfolding) fCorrelation->Fill(fill);
884 }
f2703bd2 885
5cd139bc 886 //Number of ITS cluster requirements
887 Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
888 if (recoITSnCluster){
889 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
890 icountRecoITSClusters++;
891 AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
f2703bd2 892
5cd139bc 893 Bool_t iscutsusingpid = fCuts->GetIsUsePID();
894 Int_t recoAnalysisCuts = -1, recoPidSelection = -1;
895 fCuts->SetUsePID(kFALSE);
896 recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
897
898 if (fDecayChannel==33){ // Ds case, where more possibilities are considered
899 Bool_t keepDs=ProcessDs(recoAnalysisCuts);
900 if(keepDs) recoAnalysisCuts=3;
901 }
902 else if (fDecayChannel==22){ // Lc->V0+bachelor case, where more possibilities are considered
903 Bool_t keepLctoV0bachelor=ProcessLctoV0Bachelor(recoAnalysisCuts);
904 if (keepLctoV0bachelor) recoAnalysisCuts=3;
905 }
906
6c62cb59 907
379592af 908
5cd139bc 909 fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object
910 Bool_t tempAn=(recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart);
911 if (fDecayChannel == 32) tempAn=(recoAnalysisCuts >0 || recoAnalysisCuts == isPartOrAntipart);
fdd5a897 912
5cd139bc 913 if (tempAn){
914 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight);
915 icountRecoPPR++;
916 AliDebug(3,"Reco Analysis cuts passed and container filled \n");
917 //pid selection
918 //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID);
919 //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){
920 recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
921
922 if (fDecayChannel==33){ // Ds case, where more possibilities are considered
923 Bool_t keepDs=ProcessDs(recoPidSelection);
924 if(keepDs) recoPidSelection=3;
925 } else if (fDecayChannel==22){ // Lc->V0+bachelor case, where more possibilities are considered
04b11915 926 Bool_t keepLctoV0bachelor=ProcessLctoV0Bachelor(recoPidSelection);
927 if (keepLctoV0bachelor) recoPidSelection=3;
5cd139bc 928 }
929
930 Bool_t tempPid=(recoPidSelection == 3 || recoPidSelection == isPartOrAntipart);
931 if (fDecayChannel == 32) tempPid=(recoPidSelection >0 || recoPidSelection == isPartOrAntipart);
932
933 if (tempPid){
934 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight);
935 icountRecoPID++;
936 AliDebug(3,"Reco PID cuts passed and container filled \n");
937 if(!fAcceptanceUnf){
938 Double_t fill[4]; //fill response matrix
04b11915 939 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fPDGcode,fill);
5cd139bc 940 if (bUnfolding) fCorrelation->Fill(fill);
379592af 941 }
5cd139bc 942 }
943 else {
944 AliDebug(3, "Analysis Cuts step not passed \n");
945 continue;
946 }
947 }
948 else {
949 AliDebug(3, "PID selection not passed \n");
950 continue;
951 }
952 }
953 else{
954 AliDebug(3, "Number of ITS cluster step not passed\n");
955 continue;
956 }
957 }
958 else{
959 AliDebug(3, "Reco acceptance step not passed\n");
960 continue;
961 }
962 }
963 else {
964 AliDebug(3, "Reco step not passed\n");
965 continue;
966 }
967 }
379592af 968
5cd139bc 969 if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
970 } // end loop on candidate
f2703bd2 971
5cd139bc 972 fCountReco+= icountReco;
973 fCountRecoAcc+= icountRecoAcc;
974 fCountRecoITSClusters+= icountRecoITSClusters;
975 fCountRecoPPR+= icountRecoPPR;
976 fCountRecoPID+= icountRecoPID;
379592af 977
78471b33 978 fHistEventsProcessed->Fill(1.5);
5cd139bc 979
980 delete[] containerInput;
981 delete[] containerInputMC;
982 delete cfVtxHF;
983 if (trackCuts){
984 // for (Int_t i=0; i<cfVtxHF->GetNProngs(); i++){
985 // delete [] trackCuts[i];
986 // }
987 delete [] trackCuts;
988 }
ec5288c3 989
990
379592af 991}
992
993//___________________________________________________________________________
994void AliCFTaskVertexingHF::Terminate(Option_t*)
995{
5cd139bc 996 // The Terminate() function is the last function to be called during
997 // a query. It always runs on the client, it can be used to present
998 // the results graphically or save the results to file.
379592af 999
5cd139bc 1000 AliAnalysisTaskSE::Terminate();
1001
1002 AliInfo(Form("Found %i MC particles that are %s in MC, in %d events",fCountMC,fPartName.Data(),fEvents));
1003 AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, in %d events",fCountAcc,fPartName.Data(),fEvents));
1004 AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fPartName.Data(),fEvents));
1005 AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, and satisfy ITS+TPC refit requirementin %d events",fCountRefit,fPartName.Data(),fEvents));
1006 AliInfo(Form("Found %i reco %s that are decaying in %s, in %d events",fCountReco,fPartName.Data(),fDauNames.Data(),fEvents));
1007 AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and are in the requested acceptance, in %d events",fCountRecoAcc,fPartName.Data(),fDauNames.Data(),fEvents));
1008 AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and have at least 5 clusters in ITS, in %d events",fCountRecoITSClusters,fPartName.Data(),fDauNames.Data(),fEvents));
1009 AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and satisfy PPR cuts, in %d events",fCountRecoPPR,fPartName.Data(),fDauNames.Data(),fEvents));
1010 AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and satisfy PPR+PID cuts, in %d events",fCountRecoPID,fPartName.Data(),fDauNames.Data(),fEvents));
379592af 1011
5cd139bc 1012 // draw some example plots....
1013 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
1014 if(!cont) {
1015 printf("CONTAINER NOT FOUND\n");
1016 return;
1017 }
1018 // projecting the containers to obtain histograms
1019 // first argument = variable, second argument = step
ec5288c3 1020
5cd139bc 1021 TH1D** h = new TH1D*[3];
1022 Int_t nvarToPlot = 0;
1023 if (fConfiguration == kSnail){
1024 //h = new TH1D[3][12];
1025 for (Int_t ih = 0; ih<3; ih++){
1026 if(fDecayChannel==22){
1027 nvarToPlot = 16;
1028 } else {
1029 nvarToPlot = 12;
1030 }
1031 h[ih] = new TH1D[nvarToPlot];
1032 }
1033 for(Int_t iC=1;iC<nvarToPlot; iC++){
1034 // MC-level
1035 h[0][iC] = *(cont->ShowProjection(iC,0));
1036 // MC-Acceptance level
1037 h[1][iC] = *(cont->ShowProjection(iC,1));
1038 // Reco-level
1039 h[2][iC] = *(cont->ShowProjection(iC,4));
1040 }
1041 }
1042 else {
1043 //h = new TH1D[3][12];
1044 nvarToPlot = 8;
1045 for (Int_t ih = 0; ih<3; ih++){
1046 h[ih] = new TH1D[nvarToPlot];
1047 }
1048 for(Int_t iC=0;iC<nvarToPlot; iC++){
1049 // MC-level
1050 h[0][iC] = *(cont->ShowProjection(iC,0));
1051 // MC-Acceptance level
1052 h[1][iC] = *(cont->ShowProjection(iC,1));
1053 // Reco-level
1054 h[2][iC] = *(cont->ShowProjection(iC,4));
1055 }
1056 }
1057 TString* titles;
1058 //Int_t nvarToPlot = 0;
1059 if (fConfiguration == kSnail){
1060 if(fDecayChannel==31){
1061 //nvarToPlot = 12;
1062 titles = new TString[nvarToPlot];
1063 titles[0]="pT_Dplus (GeV/c)";
1064 titles[1]="rapidity";
1065 titles[2]="phi (rad)";
1066 titles[3]="cT (#mum)";
1067 titles[4]="cosPointingAngle";
1068 titles[5]="pT_1 (GeV/c)";
1069 titles[6]="pT_2 (GeV/c)";
1070 titles[7]="pT_3 (GeV/c)";
1071 titles[8]="d0_1 (#mum)";
1072 titles[9]="d0_2 (#mum)";
1073 titles[10]="d0_3 (#mum)";
1074 titles[11]="zVertex (cm)";
1075 } else if (fDecayChannel==22) {
1076 //nvarToPlot = 16;
1077 titles = new TString[nvarToPlot];
1078 titles[0]="pT_Lc (GeV/c)";
1079 titles[1]="rapidity";
1080 titles[2]="phi (rad)";
1081 titles[3]="cosPAV0";
1082 titles[4]="onTheFlyStatusV0";
1083 titles[5]="centrality";
1084 titles[6]="fake";
1085 titles[7]="multiplicity";
1086 titles[8]="pT_bachelor (GeV/c)";
1087 titles[9]="pT_V0pos (GeV/c)";
1088 titles[10]="pT_V0neg (GeV/c)";
1089 titles[11]="invMassV0 (GeV/c2)";
1090 titles[12]="dcaV0 (nSigma)";
1091 titles[13]="c#tauV0 (#mum)";
1092 titles[14]="c#tau (#mum)";
1093 titles[15]="cosPA";
1094 } else {
1095 //nvarToPlot = 12;
1096 titles = new TString[nvarToPlot];
1097 titles[0]="pT_D0 (GeV/c)";
1098 titles[1]="rapidity";
1099 titles[2]="cosThetaStar";
1100 titles[3]="pT_pi (GeV/c)";
1101 titles[4]="pT_K (Gev/c)";
1102 titles[5]="cT (#mum)";
1103 titles[6]="dca (#mum)";
1104 titles[7]="d0_pi (#mum)";
1105 titles[8]="d0_K (#mum)";
1106 titles[9]="d0xd0 (#mum^2)";
1107 titles[10]="cosPointingAngle";
1108 titles[11]="phi (rad)";
1109 }
1110 }
1111 else {
1112 //nvarToPlot = 8;
1113 titles = new TString[nvarToPlot];
1114 if (fDecayChannel==22) {
1115 titles[0]="pT_candidate (GeV/c)";
1116 titles[1]="rapidity";
1117 titles[2]="phi (rad)";
1118 titles[3]="cosPAV0";
1119 titles[4]="onTheFlyStatusV0";
1120 titles[5]="centrality";
1121 titles[6]="fake";
1122 titles[7]="multiplicity";
1123 } else {
1124 titles[0]="pT_candidate (GeV/c)";
1125 titles[1]="rapidity";
1126 titles[2]="cT (#mum)";
1127 titles[3]="phi";
1128 titles[4]="z_{vtx}";
1129 titles[5]="centrality";
1130 titles[6]="fake";
1131 titles[7]="multiplicity";
1132 }
1133 }
1134
1135 Int_t markers[16]={20,24,21,25,27,28,
1136 20,24,21,25,27,28,
1137 20,24,21,25};
1138 Int_t colors[3]={2,8,4};
1139 for(Int_t iC=0;iC<nvarToPlot; iC++){
1140 for(Int_t iStep=0;iStep<3;iStep++){
1141 h[iStep][iC].SetTitle(titles[iC].Data());
1142 h[iStep][iC].GetXaxis()->SetTitle(titles[iC].Data());
1143 Double_t maxh=h[iStep][iC].GetMaximum();
1144 h[iStep][iC].GetYaxis()->SetRangeUser(0,maxh*1.2);
1145 h[iStep][iC].SetMarkerStyle(markers[iC]);
1146 h[iStep][iC].SetMarkerColor(colors[iStep]);
1147 }
1148 }
379592af 1149
5cd139bc 1150 gStyle->SetCanvasColor(0);
1151 gStyle->SetFrameFillColor(0);
1152 gStyle->SetTitleFillColor(0);
1153 gStyle->SetStatColor(0);
379592af 1154
5cd139bc 1155 // drawing in 2 separate canvas for a matter of clearity
1156 TCanvas * c1 =new TCanvas(Form("c1New_%d",fDecayChannel),"Vars 0, 1, 2, 3",1100,1200);
1157 c1->Divide(3,4);
1158 Int_t iPad=1;
1159 for(Int_t iVar=0; iVar<4; iVar++){
1160 c1->cd(iPad++);
1161 h[0][iVar].DrawCopy("p");
1162 c1->cd(iPad++);
1163 h[1][iVar].DrawCopy("p");
1164 c1->cd(iPad++);
1165 h[2][iVar].DrawCopy("p");
1166 }
379592af 1167
5cd139bc 1168 TCanvas * c2 =new TCanvas(Form("c2New_%d",fDecayChannel),"Vars 4, 5, 6, 7",1100,1200);
1169 c2->Divide(3,4);
1170 iPad=1;
1171 for(Int_t iVar=4; iVar<8; iVar++){
1172 c2->cd(iPad++);
1173 h[0][iVar].DrawCopy("p");
1174 c2->cd(iPad++);
1175 h[1][iVar].DrawCopy("p");
1176 c2->cd(iPad++);
1177 h[2][iVar].DrawCopy("p");
1178 }
31da6b05 1179
5cd139bc 1180 if (fConfiguration == kSnail){
1181 TCanvas * c3 =new TCanvas(Form("c3New_%d",fDecayChannel),"Vars 8, 9, 10, 11",1100,1200);
1182 c3->Divide(3,4);
1183 iPad=1;
1184 for(Int_t iVar=8; iVar<12; iVar++){
1185 c3->cd(iPad++);
1186 h[0][iVar].DrawCopy("p");
1187 c3->cd(iPad++);
1188 h[1][iVar].DrawCopy("p");
1189 c3->cd(iPad++);
1190 h[2][iVar].DrawCopy("p");
1191 }
1192 if (fDecayChannel==22) {
1193 TCanvas * c4 =new TCanvas(Form("c4New_%d",fDecayChannel),"Vars 12, 13, 14, 15",1100,1200);
1194 c4->Divide(3,4);
1195 iPad=1;
1196 for(Int_t iVar=12; iVar<16; iVar++){
1197 c4->cd(iPad++);
1198 h[0][iVar].DrawCopy("p");
1199 c4->cd(iPad++);
1200 h[1][iVar].DrawCopy("p");
1201 c4->cd(iPad++);
1202 h[2][iVar].DrawCopy("p");
1203 }
1204 }
1205 }
31da6b05 1206
5cd139bc 1207 /*
1208 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
379592af 1209
5cd139bc 1210 TH2D* corr1 = hcorr->Projection(0,2);
1211 TH2D* corr2 = hcorr->Projection(1,3);
379592af 1212
5cd139bc 1213 TCanvas * c7 =new TCanvas(Form("c7New_%d",fDecayChannel),"",800,400);
1214 c7->Divide(2,1);
1215 c7->cd(1);
1216 corr1->DrawCopy("text");
1217 c7->cd(2);
1218 corr2->DrawCopy("text");
1219 */
1220 TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
379592af 1221
5cd139bc 1222 // corr1->Write();
1223 // corr2->Write();
4e0af875 1224
5cd139bc 1225 for(Int_t iC=0;iC<nvarToPlot; iC++){
1226 for(Int_t iStep=0;iStep<3;iStep++){
1227 h[iStep][iC].Write(Form("Step%d_%s",iStep,titles[iC].Data()));
1228 }
1229 }
1230 file_projection->Close();
1231 for (Int_t ih = 0; ih<3; ih++) delete [] h[ih];
1232 delete [] h;
6299d56c 1233
379592af 1234
379592af 1235}
1236
1237//___________________________________________________________________________
f2703bd2 1238void AliCFTaskVertexingHF::UserCreateOutputObjects()
1239{
5cd139bc 1240 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1241 //TO BE SET BEFORE THE EXECUTION OF THE TASK
1242 //
1243 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
379592af 1244
5cd139bc 1245 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
1246 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1247 AliPIDResponse *localPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();
1248
1249 if (fCuts->GetIsUsePID() && fDecayChannel==22) {
1250
1251 fCuts->GetPidHF()->SetPidResponse(localPIDResponse);
5cd139bc 1252 fCuts->GetPidHF()->SetOldPid(kFALSE);
ea6f40e6 1253 AliRDHFCutsLctoV0* lcv0Cuts=dynamic_cast<AliRDHFCutsLctoV0*>(fCuts);
1254 if(lcv0Cuts){
1255 lcv0Cuts->GetPidV0pos()->SetPidResponse(localPIDResponse);
1256 lcv0Cuts->GetPidV0neg()->SetPidResponse(localPIDResponse);
1257 lcv0Cuts->GetPidV0pos()->SetOldPid(kFALSE);
1258 lcv0Cuts->GetPidV0neg()->SetOldPid(kFALSE);
1259 }
5cd139bc 1260 }
0e355b7f 1261
5cd139bc 1262 //slot #1
1263 OpenFile(1);
1264 const char* nameoutput=GetOutputSlot(1)->GetContainer()->GetName();
78471b33 1265 fHistEventsProcessed = new TH1I(nameoutput,"",2,0,2) ;
1266 fHistEventsProcessed->GetXaxis()->SetBinLabel(1,"Events processed (all)");
1267 fHistEventsProcessed->GetXaxis()->SetBinLabel(2,"Events analyzed (after selection)");
5cd139bc 1268
1269 PostData(1,fHistEventsProcessed) ;
1270 PostData(2,fCFManager->GetParticleContainer()) ;
1271 PostData(3,fCorrelation) ;
0e355b7f 1272
379592af 1273}
1274
f2703bd2 1275
942a78ce 1276//_________________________________________________________________________
1277void AliCFTaskVertexingHF::SetPtWeightsFromFONLL276overLHC12a17a(){
1278 // ad-hoc weight function from ratio of
1279 // pt spectra from FONLL 2.76 TeV and
1280 // pt spectra from MC production LHC12a17a (PYTHIA Perugia0 with pthard bins)
1281 if(fFuncWeight) delete fFuncWeight;
1282 fFuncWeight=new TF1("funcWeight","[0]+[1]*TMath::Exp(-[2]*x)",0.,50.);
1283 fFuncWeight->SetParameter(0,4.63891e-02);
1284 fFuncWeight->SetParameter(1,1.51674e+01);
1285 fFuncWeight->SetParameter(2,4.09941e-01);
1286}
1287//_________________________________________________________________________
1288void AliCFTaskVertexingHF::SetPtWeightsFromDataPbPb276overLHC12a17a(){
1289 // ad-hoc weight function from ratio of
1290 // D0 pt spectra in PbPb 2011 0-10% centrality and
1291 // pt spectra from MC production LHC12a17a (PYTHIA Perugia0 with pthard bins)
1292 if(fFuncWeight) delete fFuncWeight;
1293 fFuncWeight=new TF1("funcWeight","[0]+[1]/TMath::Power(x,[2])",0.05,50.);
1294 fFuncWeight->SetParameter(0,1.43116e-02);
1295 fFuncWeight->SetParameter(1,4.37758e+02);
1296 fFuncWeight->SetParameter(2,3.08583);
1297}
f2703bd2 1298//_________________________________________________________________________
1299Double_t AliCFTaskVertexingHF::GetWeight(Float_t pt)
1300{
5cd139bc 1301 //
1302 // calculating the weight to fill the container
1303 //
f2703bd2 1304
5cd139bc 1305 // FNOLL central:
1306 // p0 = 1.63297e-01 --> 0.322643
1307 // p1 = 2.96275e+00
1308 // p2 = 2.30301e+00
1309 // p3 = 2.50000e+00
1310
1311 // PYTHIA
1312 // p0 = 1.85906e-01 --> 0.36609
1313 // p1 = 1.94635e+00
1314 // p2 = 1.40463e+00
1315 // p3 = 2.50000e+00
1316
1317 Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
1318 Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
1319
1320 Double_t dndpt_func1 = dNdptFit(pt,func1);
1321 if(fUseFlatPtWeight) dndpt_func1 = 1./30.;
1322 Double_t dndpt_func2 = dNdptFit(pt,func2);
1323 AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
1324 return dndpt_func1/dndpt_func2;
f2703bd2 1325}
1326
1327//__________________________________________________________________________________________________
1328Double_t AliCFTaskVertexingHF::dNdptFit(Float_t pt, Double_t* par)
1329{
5cd139bc 1330 //
1331 // calculating dNdpt
1332 //
f2703bd2 1333
5cd139bc 1334 Double_t denom = TMath::Power((pt/par[1]), par[3] );
1335 Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
f2703bd2 1336
5cd139bc 1337 return dNdpt;
f2703bd2 1338}
6c62cb59 1339
17bf1a34 1340//__________________________________________________________________________________________________
1341Double_t AliCFTaskVertexingHF::GetZWeight(Float_t z, Int_t runnumber){
1342 //
1343 // calculating the z-vtx weight for the given run range
1344 //
1345
1346 if(runnumber>146824 || runnumber<146803) return 1.0;
1347
1348 Double_t func1[3] = {1.0, -0.5, 6.5 };
1349 Double_t func2[3] = {1.0, -0.5, 5.5 };
1350
1351 Double_t dzFunc1 = DodzFit(z,func1);
1352 Double_t dzFunc2 = DodzFit(z,func2);
1353
1354 return dzFunc1/dzFunc2;
1355}
1356
1357//__________________________________________________________________________________________________
1358Double_t AliCFTaskVertexingHF::DodzFit(Float_t z, Double_t* par) {
1359
1360 //
1361 // Gaussian z-vtx shape
1362 //
1363 //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
1364
1365 Double_t value = par[0]/TMath::Sqrt(2.*TMath::Pi())/par[2]*TMath::Exp(-(z-par[1])*(z-par[1])/2./par[2]/par[2]);
1366
1367 return value;
1368}
6c62cb59 1369//__________________________________________________________________________________________________
afb24c40 1370Double_t AliCFTaskVertexingHF::GetNchWeight(Int_t nch){
1371 //
1372 // calculates the Nch weight using the measured and generateed Nch distributions
1373 //
1374 if(nch<=0) return 0.;
1375 Double_t pMeas=fHistoMeasNch->GetBinContent(fHistoMeasNch->FindBin(nch));
1376 Double_t pMC=fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(nch));
1377 return pMeas/pMC;
1378}
1379//__________________________________________________________________________________________________
1380void AliCFTaskVertexingHF::CreateMeasuredNchHisto(){
1381 // creates historgam with measured multiplcity distribution in pp 7 TeV collisions (from Eur. Phys. J. C (2010) 68: 345–354)
1382 Double_t nchbins[66]={0.50,1.50,2.50,3.50,4.50,5.50,6.50,7.50,8.50,9.50,
1383 10.50,11.50,12.50,13.50,14.50,15.50,16.50,17.50,18.50,19.50,
1384 20.50,21.50,22.50,23.50,24.50,25.50,26.50,27.50,28.50,29.50,
1385 30.50,31.50,32.50,33.50,34.50,35.50,36.50,37.50,38.50,39.50,
1386 40.50,41.50,42.50,43.50,44.50,45.50,46.50,47.50,48.50,49.50,
1387 50.50,51.50,52.50,53.50,54.50,55.50,56.50,57.50,58.50,59.50,
1388 60.50,62.50,64.50,66.50,68.50,70.50};
1389 Double_t pch[65]={0.062011,0.072943,0.070771,0.067245,0.062834,0.057383,0.051499,0.04591,0.041109,0.036954,
1390 0.03359,0.030729,0.028539,0.026575,0.024653,0.0229,0.021325,0.019768,0.018561,0.017187,
1391 0.01604,0.014836,0.013726,0.012576,0.011481,0.010393,0.009502,0.008776,0.008024,0.007452,
1392 0.006851,0.006428,0.00594,0.005515,0.005102,0.00469,0.004162,0.003811,0.003389,0.003071,
1393 0.002708,0.002422,0.002184,0.001968,0.00186,0.00165,0.001577,0.001387,0.001254,0.001118,
1394 0.001037,0.000942,0.000823,0.000736,0.000654,0.000579,0.000512,0.00049,0.00045,0.000355,
1395 0.000296,0.000265,0.000193,0.00016,0.000126};
1396
1397 if(fHistoMeasNch) delete fHistoMeasNch;
1398 fHistoMeasNch=new TH1F("hMeaseNch","",65,nchbins);
1399 for(Int_t i=0; i<65; i++){
1400 fHistoMeasNch->SetBinContent(i+1,pch[i]);
1401 fHistoMeasNch->SetBinError(i+1,0.);
1402 }
1403}
1404
5cd139bc 1405//__________________________________________________________________________________________________
6c62cb59 1406Bool_t AliCFTaskVertexingHF::ProcessDs(Int_t recoAnalysisCuts) const{
1407 // processes output of Ds is selected
1408 Bool_t keep=kFALSE;
1409 if(recoAnalysisCuts > 0){
1410 Int_t isKKpi=recoAnalysisCuts&1;
1411 Int_t ispiKK=recoAnalysisCuts&2;
1412 Int_t isPhiKKpi=recoAnalysisCuts&4;
1413 Int_t isPhipiKK=recoAnalysisCuts&8;
1414 Int_t isK0starKKpi=recoAnalysisCuts&16;
1415 Int_t isK0starpiKK=recoAnalysisCuts&32;
1416 if(fDsOption==1){
1417 if(isKKpi && isPhiKKpi) keep=kTRUE;
1418 if(ispiKK && isPhipiKK) keep=kTRUE;
1419 }
1420 else if(fDsOption==2){
1421 if(isKKpi && isK0starKKpi) keep=kTRUE;
1422 if(ispiKK && isK0starpiKK) keep=kTRUE;
1423 }
1424 else if(fDsOption==3)keep=kTRUE;
1425 }
1426 return keep;
1427}
5cd139bc 1428//__________________________________________________________________________________________________
1429Bool_t AliCFTaskVertexingHF::ProcessLctoV0Bachelor(Int_t recoAnalysisCuts) const{
1430 // processes output of Lc->V0+bnachelor is selected
1431 Bool_t keep=kFALSE;
1432 if(recoAnalysisCuts > 0){
1433
1434 Int_t isK0Sp = recoAnalysisCuts&1;
1435 Int_t isLambdaBarpi = recoAnalysisCuts&2;
1436 Int_t isLambdapi = recoAnalysisCuts&4;
1437
1438 if(fLctoV0bachelorOption==1){
1439 if(isK0Sp) keep=kTRUE;
1440 }
1441 else if(fLctoV0bachelorOption==2){
1442 if(isLambdaBarpi) keep=kTRUE;
1443 }
1444 else if(fLctoV0bachelorOption==4){
1445 if(isLambdapi) keep=kTRUE;
1446 }
1447 else if(fLctoV0bachelorOption==7) {
1448 if (isK0Sp || isLambdaBarpi || isLambdapi) keep=kTRUE;
1449 }
1450 }
1451 return keep;
1452}