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