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