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