]> 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
from Andrea to check 10-20% v2
[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;
686 return;
687 }
751b8274 688
8bbb4282 689 if(aodEvent->GetTriggerMask()==0 &&
690 (runnumber>=195344 && runnumber<=195677)){
691 AliDebug(3,"Event rejected because of null trigger mask");
692 delete[] containerInput;
693 delete[] containerInputMC;
694 return;
695 }
696
5cd139bc 697 AliESDtrackCuts** trackCuts = new AliESDtrackCuts*[cfVtxHF->GetNProngs()];
698 if (fDecayChannel == 21){
699 // for the D*, setting the third element of the array of the track cuts to those for the soft pion
700 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs()-1; iProng++){
701 trackCuts[iProng]=fCuts->GetTrackCuts();
702 }
703 trackCuts[2] = fCuts->GetTrackCutsSoftPi();
704 }
705 else if (fDecayChannel == 22) {
706 // for the Lc->V0+bachelor, setting the second and third elements of the array of the track cuts to those for the V0 daughters
707 trackCuts[0]=fCuts->GetTrackCuts();
708 trackCuts[1]=fCuts->GetTrackCutsV0daughters();
709 trackCuts[2]=fCuts->GetTrackCutsV0daughters();
710 }
711 else {
712 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs(); iProng++){
713 trackCuts[iProng]=fCuts->GetTrackCuts();
714 }
715 }
2bf2e62b 716
5cd139bc 717 //General settings: vertex, feed down and fill reco container with generated values.
718 cfVtxHF->SetRecoPrimVertex(zPrimVertex);
719 cfVtxHF->SetMCPrimaryVertex(zMCVertex);
720 cfVtxHF->SetFillFromGenerated(fFillFromGenerated);
721 cfVtxHF->SetNVar(fNvar);
722 cfVtxHF->SetFakeSelection(fFakeSelection);
723 cfVtxHF->SetRejectCandidateIfNotFromQuark(fRejectIfNoQuark);
724 cfVtxHF->SetConfiguration(fConfiguration);
725
726 // switch-off the trigger class selection (doesn't work for MC)
727 fCuts->SetTriggerClass("");
728
729 // MC vertex, to be used, in case, for pp
730 if (fUseMCVertex) fCuts->SetUseMCVertex();
731
732 if (fCentralitySelection){ // keep only the requested centrality
733 if(fCuts->IsEventSelectedInCentrality(aodEvent)!=0) {
734 delete[] containerInput;
735 delete[] containerInputMC;
736 delete [] trackCuts;
737 return;
738 }
739 } else { // keep all centralities
740 fCuts->SetMinCentrality(0.);
741 fCuts->SetMaxCentrality(100.);
742 }
8b1b055d 743
d74f7bc1 744 Float_t centValue = 0.;
745 if(!fIsPPData) centValue = fCuts->GetCentrality(aodEvent);
5cd139bc 746 cfVtxHF->SetCentralityValue(centValue);
f2703bd2 747
f26726ed 748 // multiplicity estimator with VZERO
d74f7bc1 749 Double_t vzeroMult=0;
750 AliAODVZERO *vzeroAOD = (AliAODVZERO*)aodEvent->GetVZEROData();
751 if(vzeroAOD) vzeroMult = vzeroAOD->GetMTotV0A() + vzeroAOD->GetMTotV0C();
752
f26726ed 753 Double_t multiplicity = nTracklets; // set to the Ntracklet estimator
d74f7bc1 754 if(fMultiplicityEstimator==kVZERO) { multiplicity = vzeroMult; }
755
756
5cd139bc 757 cfVtxHF->SetMultiplicity(multiplicity);
d74f7bc1 758
759 // printf("Multiplicity estimator %d, value %2.2f\n",fMultiplicityEstimator,multiplicity);
ec5288c3 760
5cd139bc 761 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
762 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
763 if (!mcPart){
764 AliError("Failed casting particle from MC array!, Skipping particle");
765 continue;
766 }
04b11915 767
768 //counting c quarks
769 cquarks += cfVtxHF->MCcquarkCounting(mcPart);
770
5cd139bc 771 // check the MC-level cuts, must be the desidered particle
772 if (!fCFManager->CheckParticleCuts(0, mcPart)) {
773 AliDebug(2,"Check the MC-level cuts - not desidered particle");
774 continue; // 0 stands for MC level
775 }
776 cfVtxHF->SetMCCandidateParam(iPart);
b7af2639 777
b7af2639 778
5cd139bc 779 if (!(cfVtxHF->SetLabelArray())){
780 AliDebug(2,Form("Impossible to set the label array (decaychannel = %d)",fDecayChannel));
781 continue;
782 }
2bf2e62b 783
5cd139bc 784 //check the candiate family at MC level
785 if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) {
786 AliDebug(2,Form("Check on the family wrong!!! (decaychannel = %d)",fDecayChannel));
787 continue;
788 }
789 else{
d74f7bc1 790 AliDebug(2,Form("Check on the family OK!!! (decaychannel = %d)",fDecayChannel));
5cd139bc 791 }
f2703bd2 792
5cd139bc 793 //Fill the MC container
794 Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
795 AliDebug(2,Form("mcContainerFilled = %d)",mcContainerFilled));
796 if (mcContainerFilled) {
797 if (fUseWeight){
798 if (fFuncWeight){ // using user-defined function
799 AliDebug(2,"Using function");
f26726ed 800 fWeight = eventWeight*fFuncWeight->Eval(containerInputMC[0]);
5cd139bc 801 }
802 else{ // using FONLL
803 AliDebug(2,"Using FONLL");
f26726ed 804 fWeight = eventWeight*GetWeight(containerInputMC[0]);
5cd139bc 805 }
806 AliDebug(2,Form("pt = %f, weight = %f",containerInputMC[0], fWeight));
807 }
808 if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) continue;
809 //MC Limited Acceptance
810 if (TMath::Abs(containerInputMC[1]) < 0.5) {
811 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight);
812 AliDebug(3,"MC Lim Acc container filled\n");
813 }
f2703bd2 814
5cd139bc 815 //MC
816 fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight);
817 icountMC++;
818 AliDebug(3,"MC cointainer filled \n");
f2703bd2 819
5cd139bc 820 // MC in acceptance
821 // check the MC-Acceptance level cuts
822 // since standard CF functions are not applicable, using Kine Cuts on daughters
823 Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep();
824 if (mcAccepStep){
825 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight);
826 AliDebug(3,"MC acceptance cut passed\n");
827 icountAcc++;
f2703bd2 828
5cd139bc 829 //MC Vertex step
830 if (fCuts->IsEventSelected(aodEvent)){
831 // filling the container if the vertex is ok
832 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex, fWeight) ;
833 AliDebug(3,"Vertex cut passed and container filled\n");
834 icountVertex++;
f2703bd2 835
5cd139bc 836 //mc Refit requirement
837 Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, &trackCuts[0]);
838 if (mcRefitStep){
839 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight);
840 AliDebug(3,"MC Refit cut passed and container filled\n");
841 icountRefit++;
842 }
843 else{
844 AliDebug(3,"MC Refit cut not passed\n");
845 continue;
846 }
379592af 847 }
5cd139bc 848 else{
849 AliDebug (3, "MC vertex step not passed\n");
850 continue;
851 }
852 }
853 else{
854 AliDebug (3, "MC in acceptance step not passed\n");
855 continue;
856 }
857 }
858 else {
859 AliDebug (3, "MC container not filled\n");
860 }
861 }
f2703bd2 862
5cd139bc 863 if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks));
864 AliDebug(2,Form("Found %i MC particles that are %s!!",icountMC,fPartName.Data()));
865 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Acc cuts!!",icountAcc,fPartName.Data()));
866 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Vertex cuts!!",icountVertex,fPartName.Data()));
867 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Refit cuts!!",icountRefit,fPartName.Data()));
868
869 // Now go to rec level
870 fCountMC += icountMC;
871 fCountAcc += icountAcc;
872 fCountVertex+= icountVertex;
873 fCountRefit+= icountRefit;
874
875 AliDebug(2,Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel));
379592af 876
5cd139bc 877 for(Int_t iCandid = 0; iCandid<arrayBranch->GetEntriesFast();iCandid++){
878 AliAODRecoDecayHF* charmCandidate=0x0;
879 switch (fDecayChannel){
880 case 2:{
881 charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid);
882 break;
883 }
884 case 21:
885 case 22:{
886 charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid);
887 break;
888 }
889 case 31:
890 case 32:
891 case 33:{
892 charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid);
893 break;
894 }
895 case 4:{
896 charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid);
897 break;
898 }
899 default:
900 break;
901 }
f2703bd2 902
5cd139bc 903 Bool_t unsetvtx=kFALSE;
904 if(!charmCandidate->GetOwnPrimaryVtx()) {
905 charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
906 unsetvtx=kTRUE;
907 }
f2703bd2 908
5cd139bc 909 Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate);
910 if (!signAssociation){
911 charmCandidate = 0x0;
912 continue;
913 }
3ee5eb83 914
5cd139bc 915 Int_t isPartOrAntipart = cfVtxHF->CheckReflexion(fSign);
916 if (isPartOrAntipart == 0){
917 AliDebug(2, Form("The candidate pdg code doesn't match the requirement set in the task (fSign = %d)",fSign));
918 continue;
919 }
3ee5eb83 920
5cd139bc 921 AliDebug(3,Form("iCandid=%d - signAssociation=%d, isPartOrAntipart=%d",iCandid, signAssociation, isPartOrAntipart));
6694a2a8 922
5cd139bc 923 Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
924 if (recoContFilled){
6694a2a8 925
5cd139bc 926 // weight according to pt
927 if (fUseWeight){
928 if (fFuncWeight){ // using user-defined function
929 AliDebug(2, "Using function");
f26726ed 930 fWeight = eventWeight*fFuncWeight->Eval(containerInput[0]);
5cd139bc 931 }
932 else{ // using FONLL
933 AliDebug(2, "Using FONLL");
f26726ed 934 fWeight = eventWeight*GetWeight(containerInput[0]);
5cd139bc 935 }
936 AliDebug(2, Form("pt = %f, weight = %f",containerInput[0], fWeight));
937 }
6694a2a8 938
5cd139bc 939 if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue;
f2703bd2 940
5cd139bc 941 //Reco Step
942 Bool_t recoStep = cfVtxHF->RecoStep();
943 Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
f2703bd2 944
5cd139bc 945
ca7d039b 946 // Selection on the filtering bit
947 Bool_t isBitSelected = kTRUE;
948 if(fDecayChannel==2) {
949 if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)) isBitSelected = kFALSE;
950 }else if(fDecayChannel==31){
951 if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kDplusCuts)) isBitSelected = kFALSE;
952 }else if(fDecayChannel==33){
953 if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kDsCuts)) isBitSelected = kFALSE;
954 }
955 if(!isBitSelected) continue;
956
957
958
5cd139bc 959 if (recoStep && recoContFilled && vtxCheck){
960 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
961 icountReco++;
962 AliDebug(3,"Reco step passed and container filled\n");
f2703bd2 963
5cd139bc 964 //Reco in the acceptance -- take care of UNFOLDING!!!!
965 Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(&trackCuts[0]);
966 if (recoAcceptanceStep) {
967 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
968 icountRecoAcc++;
969 AliDebug(3,"Reco acceptance cut passed and container filled\n");
f2703bd2 970
5cd139bc 971 if(fAcceptanceUnf){
972 Double_t fill[4]; //fill response matrix
04b11915 973 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fPDGcode,fill);
5cd139bc 974 if (bUnfolding) fCorrelation->Fill(fill);
975 }
f2703bd2 976
5cd139bc 977 //Number of ITS cluster requirements
978 Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
979 if (recoITSnCluster){
980 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
981 icountRecoITSClusters++;
982 AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
f2703bd2 983
5cd139bc 984 Bool_t iscutsusingpid = fCuts->GetIsUsePID();
985 Int_t recoAnalysisCuts = -1, recoPidSelection = -1;
986 fCuts->SetUsePID(kFALSE);
987 recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
988
989 if (fDecayChannel==33){ // Ds case, where more possibilities are considered
990 Bool_t keepDs=ProcessDs(recoAnalysisCuts);
991 if(keepDs) recoAnalysisCuts=3;
992 }
993 else if (fDecayChannel==22){ // Lc->V0+bachelor case, where more possibilities are considered
994 Bool_t keepLctoV0bachelor=ProcessLctoV0Bachelor(recoAnalysisCuts);
995 if (keepLctoV0bachelor) recoAnalysisCuts=3;
996 }
997
6c62cb59 998
379592af 999
5cd139bc 1000 fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object
1001 Bool_t tempAn=(recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart);
1002 if (fDecayChannel == 32) tempAn=(recoAnalysisCuts >0 || recoAnalysisCuts == isPartOrAntipart);
fdd5a897 1003
5cd139bc 1004 if (tempAn){
1005 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight);
1006 icountRecoPPR++;
1007 AliDebug(3,"Reco Analysis cuts passed and container filled \n");
1008 //pid selection
1009 //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID);
1010 //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){
1011 recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
1012
1013 if (fDecayChannel==33){ // Ds case, where more possibilities are considered
1014 Bool_t keepDs=ProcessDs(recoPidSelection);
1015 if(keepDs) recoPidSelection=3;
1016 } else if (fDecayChannel==22){ // Lc->V0+bachelor case, where more possibilities are considered
04b11915 1017 Bool_t keepLctoV0bachelor=ProcessLctoV0Bachelor(recoPidSelection);
1018 if (keepLctoV0bachelor) recoPidSelection=3;
5cd139bc 1019 }
1020
1021 Bool_t tempPid=(recoPidSelection == 3 || recoPidSelection == isPartOrAntipart);
1022 if (fDecayChannel == 32) tempPid=(recoPidSelection >0 || recoPidSelection == isPartOrAntipart);
1023
1024 if (tempPid){
1025 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight);
1026 icountRecoPID++;
1027 AliDebug(3,"Reco PID cuts passed and container filled \n");
1028 if(!fAcceptanceUnf){
1029 Double_t fill[4]; //fill response matrix
04b11915 1030 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fPDGcode,fill);
5cd139bc 1031 if (bUnfolding) fCorrelation->Fill(fill);
379592af 1032 }
5cd139bc 1033 }
1034 else {
1035 AliDebug(3, "Analysis Cuts step not passed \n");
1036 continue;
1037 }
1038 }
1039 else {
1040 AliDebug(3, "PID selection not passed \n");
1041 continue;
1042 }
1043 }
1044 else{
1045 AliDebug(3, "Number of ITS cluster step not passed\n");
1046 continue;
1047 }
1048 }
1049 else{
1050 AliDebug(3, "Reco acceptance step not passed\n");
1051 continue;
1052 }
1053 }
1054 else {
1055 AliDebug(3, "Reco step not passed\n");
1056 continue;
1057 }
1058 }
379592af 1059
5cd139bc 1060 if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1061 } // end loop on candidate
f2703bd2 1062
5cd139bc 1063 fCountReco+= icountReco;
1064 fCountRecoAcc+= icountRecoAcc;
1065 fCountRecoITSClusters+= icountRecoITSClusters;
1066 fCountRecoPPR+= icountRecoPPR;
1067 fCountRecoPID+= icountRecoPID;
379592af 1068
78471b33 1069 fHistEventsProcessed->Fill(1.5);
5cd139bc 1070
1071 delete[] containerInput;
1072 delete[] containerInputMC;
1073 delete cfVtxHF;
1074 if (trackCuts){
1075 // for (Int_t i=0; i<cfVtxHF->GetNProngs(); i++){
1076 // delete [] trackCuts[i];
1077 // }
1078 delete [] trackCuts;
1079 }
ec5288c3 1080
1081
379592af 1082}
1083
1084//___________________________________________________________________________
1085void AliCFTaskVertexingHF::Terminate(Option_t*)
1086{
5cd139bc 1087 // The Terminate() function is the last function to be called during
1088 // a query. It always runs on the client, it can be used to present
1089 // the results graphically or save the results to file.
379592af 1090
5cd139bc 1091 AliAnalysisTaskSE::Terminate();
1092
1093 AliInfo(Form("Found %i MC particles that are %s in MC, in %d events",fCountMC,fPartName.Data(),fEvents));
1094 AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, in %d events",fCountAcc,fPartName.Data(),fEvents));
1095 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));
1096 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));
1097 AliInfo(Form("Found %i reco %s that are decaying in %s, in %d events",fCountReco,fPartName.Data(),fDauNames.Data(),fEvents));
1098 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));
1099 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));
1100 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));
1101 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 1102
5cd139bc 1103 // draw some example plots....
1104 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
1105 if(!cont) {
1106 printf("CONTAINER NOT FOUND\n");
1107 return;
1108 }
1109 // projecting the containers to obtain histograms
1110 // first argument = variable, second argument = step
ec5288c3 1111
5cd139bc 1112 TH1D** h = new TH1D*[3];
1113 Int_t nvarToPlot = 0;
1114 if (fConfiguration == kSnail){
1115 //h = new TH1D[3][12];
1116 for (Int_t ih = 0; ih<3; ih++){
1117 if(fDecayChannel==22){
1118 nvarToPlot = 16;
1119 } else {
1120 nvarToPlot = 12;
1121 }
1122 h[ih] = new TH1D[nvarToPlot];
1123 }
1124 for(Int_t iC=1;iC<nvarToPlot; iC++){
1125 // MC-level
1126 h[0][iC] = *(cont->ShowProjection(iC,0));
1127 // MC-Acceptance level
1128 h[1][iC] = *(cont->ShowProjection(iC,1));
1129 // Reco-level
1130 h[2][iC] = *(cont->ShowProjection(iC,4));
1131 }
1132 }
1133 else {
1134 //h = new TH1D[3][12];
1135 nvarToPlot = 8;
1136 for (Int_t ih = 0; ih<3; ih++){
1137 h[ih] = new TH1D[nvarToPlot];
1138 }
1139 for(Int_t iC=0;iC<nvarToPlot; iC++){
1140 // MC-level
1141 h[0][iC] = *(cont->ShowProjection(iC,0));
1142 // MC-Acceptance level
1143 h[1][iC] = *(cont->ShowProjection(iC,1));
1144 // Reco-level
1145 h[2][iC] = *(cont->ShowProjection(iC,4));
1146 }
1147 }
1148 TString* titles;
1149 //Int_t nvarToPlot = 0;
1150 if (fConfiguration == kSnail){
1151 if(fDecayChannel==31){
1152 //nvarToPlot = 12;
1153 titles = new TString[nvarToPlot];
1154 titles[0]="pT_Dplus (GeV/c)";
1155 titles[1]="rapidity";
1156 titles[2]="phi (rad)";
1157 titles[3]="cT (#mum)";
1158 titles[4]="cosPointingAngle";
1159 titles[5]="pT_1 (GeV/c)";
1160 titles[6]="pT_2 (GeV/c)";
1161 titles[7]="pT_3 (GeV/c)";
1162 titles[8]="d0_1 (#mum)";
1163 titles[9]="d0_2 (#mum)";
1164 titles[10]="d0_3 (#mum)";
1165 titles[11]="zVertex (cm)";
1166 } else if (fDecayChannel==22) {
1167 //nvarToPlot = 16;
1168 titles = new TString[nvarToPlot];
53c91ed7 1169 titles[0]="p_{T}(#Lambda_{c}) [GeV/c]";
1170 titles[1]="y(#Lambda_{c})";
1171 titles[2]="#varphi(#Lambda_{c}) [rad]";
1172 titles[3]="onTheFlyStatusV0";
1173 titles[4]="z_{vtx} [cm]";
5cd139bc 1174 titles[5]="centrality";
1175 titles[6]="fake";
1176 titles[7]="multiplicity";
53c91ed7 1177 //titles[8]="pT(bachelor) [GeV/c]";
1178 titles[8]="p(bachelor) [GeV/c]";
1179 titles[9]="p_{T}(V0) [GeV/c]";
1180 titles[10]="y(V0)";
1181 titles[11]="#varphi(V0) [rad]";
1182 titles[12]="m_{inv}(#pi^{+}#pi^{+}) [GeV/c^{2}]";
1183 titles[13]="dcaV0 (nSigma)";
1184 titles[14]="cosine pointing angle (V0)";
1185 titles[15]="cosine pointing angle (#Lambda_{c})";
1186 //titles[16]="c#tauV0 (#mum)";
1187 //titles[17]="c#tau (#mum)";
5cd139bc 1188 } else {
1189 //nvarToPlot = 12;
1190 titles = new TString[nvarToPlot];
1191 titles[0]="pT_D0 (GeV/c)";
1192 titles[1]="rapidity";
1193 titles[2]="cosThetaStar";
1194 titles[3]="pT_pi (GeV/c)";
1195 titles[4]="pT_K (Gev/c)";
1196 titles[5]="cT (#mum)";
1197 titles[6]="dca (#mum)";
1198 titles[7]="d0_pi (#mum)";
1199 titles[8]="d0_K (#mum)";
1200 titles[9]="d0xd0 (#mum^2)";
1201 titles[10]="cosPointingAngle";
1202 titles[11]="phi (rad)";
1203 }
1204 }
1205 else {
1206 //nvarToPlot = 8;
1207 titles = new TString[nvarToPlot];
1208 if (fDecayChannel==22) {
53c91ed7 1209 titles[0]="p_{T}(#Lambda_{c}) [GeV/c]";
1210 titles[1]="y(#Lambda_{c})";
1211 titles[2]="#varphi(#Lambda_{c}) [rad]";
1212 titles[3]="onTheFlyStatusV0";
1213 titles[4]="z_{vtx} [cm]";
5cd139bc 1214 titles[5]="centrality";
1215 titles[6]="fake";
1216 titles[7]="multiplicity";
1217 } else {
1218 titles[0]="pT_candidate (GeV/c)";
1219 titles[1]="rapidity";
1220 titles[2]="cT (#mum)";
1221 titles[3]="phi";
1222 titles[4]="z_{vtx}";
1223 titles[5]="centrality";
1224 titles[6]="fake";
1225 titles[7]="multiplicity";
1226 }
1227 }
1228
1229 Int_t markers[16]={20,24,21,25,27,28,
1230 20,24,21,25,27,28,
1231 20,24,21,25};
1232 Int_t colors[3]={2,8,4};
1233 for(Int_t iC=0;iC<nvarToPlot; iC++){
1234 for(Int_t iStep=0;iStep<3;iStep++){
1235 h[iStep][iC].SetTitle(titles[iC].Data());
1236 h[iStep][iC].GetXaxis()->SetTitle(titles[iC].Data());
1237 Double_t maxh=h[iStep][iC].GetMaximum();
1238 h[iStep][iC].GetYaxis()->SetRangeUser(0,maxh*1.2);
1239 h[iStep][iC].SetMarkerStyle(markers[iC]);
1240 h[iStep][iC].SetMarkerColor(colors[iStep]);
1241 }
1242 }
379592af 1243
5cd139bc 1244 gStyle->SetCanvasColor(0);
1245 gStyle->SetFrameFillColor(0);
1246 gStyle->SetTitleFillColor(0);
1247 gStyle->SetStatColor(0);
379592af 1248
5cd139bc 1249 // drawing in 2 separate canvas for a matter of clearity
1250 TCanvas * c1 =new TCanvas(Form("c1New_%d",fDecayChannel),"Vars 0, 1, 2, 3",1100,1200);
1251 c1->Divide(3,4);
1252 Int_t iPad=1;
1253 for(Int_t iVar=0; iVar<4; iVar++){
1254 c1->cd(iPad++);
1255 h[0][iVar].DrawCopy("p");
1256 c1->cd(iPad++);
1257 h[1][iVar].DrawCopy("p");
1258 c1->cd(iPad++);
1259 h[2][iVar].DrawCopy("p");
1260 }
379592af 1261
5cd139bc 1262 TCanvas * c2 =new TCanvas(Form("c2New_%d",fDecayChannel),"Vars 4, 5, 6, 7",1100,1200);
1263 c2->Divide(3,4);
1264 iPad=1;
1265 for(Int_t iVar=4; iVar<8; iVar++){
1266 c2->cd(iPad++);
1267 h[0][iVar].DrawCopy("p");
1268 c2->cd(iPad++);
1269 h[1][iVar].DrawCopy("p");
1270 c2->cd(iPad++);
1271 h[2][iVar].DrawCopy("p");
1272 }
31da6b05 1273
5cd139bc 1274 if (fConfiguration == kSnail){
1275 TCanvas * c3 =new TCanvas(Form("c3New_%d",fDecayChannel),"Vars 8, 9, 10, 11",1100,1200);
1276 c3->Divide(3,4);
1277 iPad=1;
1278 for(Int_t iVar=8; iVar<12; iVar++){
1279 c3->cd(iPad++);
1280 h[0][iVar].DrawCopy("p");
1281 c3->cd(iPad++);
1282 h[1][iVar].DrawCopy("p");
1283 c3->cd(iPad++);
1284 h[2][iVar].DrawCopy("p");
1285 }
1286 if (fDecayChannel==22) {
1287 TCanvas * c4 =new TCanvas(Form("c4New_%d",fDecayChannel),"Vars 12, 13, 14, 15",1100,1200);
1288 c4->Divide(3,4);
1289 iPad=1;
1290 for(Int_t iVar=12; iVar<16; iVar++){
1291 c4->cd(iPad++);
1292 h[0][iVar].DrawCopy("p");
1293 c4->cd(iPad++);
1294 h[1][iVar].DrawCopy("p");
1295 c4->cd(iPad++);
1296 h[2][iVar].DrawCopy("p");
1297 }
1298 }
1299 }
31da6b05 1300
5cd139bc 1301 /*
1302 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
379592af 1303
5cd139bc 1304 TH2D* corr1 = hcorr->Projection(0,2);
1305 TH2D* corr2 = hcorr->Projection(1,3);
379592af 1306
5cd139bc 1307 TCanvas * c7 =new TCanvas(Form("c7New_%d",fDecayChannel),"",800,400);
1308 c7->Divide(2,1);
1309 c7->cd(1);
1310 corr1->DrawCopy("text");
1311 c7->cd(2);
1312 corr2->DrawCopy("text");
1313 */
1314 TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
379592af 1315
5cd139bc 1316 // corr1->Write();
1317 // corr2->Write();
4e0af875 1318
5cd139bc 1319 for(Int_t iC=0;iC<nvarToPlot; iC++){
1320 for(Int_t iStep=0;iStep<3;iStep++){
1321 h[iStep][iC].Write(Form("Step%d_%s",iStep,titles[iC].Data()));
1322 }
1323 }
1324 file_projection->Close();
1325 for (Int_t ih = 0; ih<3; ih++) delete [] h[ih];
1326 delete [] h;
6299d56c 1327
379592af 1328
379592af 1329}
1330
1331//___________________________________________________________________________
f2703bd2 1332void AliCFTaskVertexingHF::UserCreateOutputObjects()
1333{
5cd139bc 1334 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1335 //TO BE SET BEFORE THE EXECUTION OF THE TASK
1336 //
1337 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
0e355b7f 1338
5cd139bc 1339 //slot #1
1340 OpenFile(1);
1341 const char* nameoutput=GetOutputSlot(1)->GetContainer()->GetName();
78471b33 1342 fHistEventsProcessed = new TH1I(nameoutput,"",2,0,2) ;
1343 fHistEventsProcessed->GetXaxis()->SetBinLabel(1,"Events processed (all)");
1344 fHistEventsProcessed->GetXaxis()->SetBinLabel(2,"Events analyzed (after selection)");
5cd139bc 1345
1346 PostData(1,fHistEventsProcessed) ;
1347 PostData(2,fCFManager->GetParticleContainer()) ;
1348 PostData(3,fCorrelation) ;
0e355b7f 1349
379592af 1350}
1351
f2703bd2 1352
942a78ce 1353//_________________________________________________________________________
1354void AliCFTaskVertexingHF::SetPtWeightsFromFONLL276overLHC12a17a(){
1355 // ad-hoc weight function from ratio of
1356 // pt spectra from FONLL 2.76 TeV and
1357 // pt spectra from MC production LHC12a17a (PYTHIA Perugia0 with pthard bins)
1358 if(fFuncWeight) delete fFuncWeight;
1359 fFuncWeight=new TF1("funcWeight","[0]+[1]*TMath::Exp(-[2]*x)",0.,50.);
1360 fFuncWeight->SetParameter(0,4.63891e-02);
1361 fFuncWeight->SetParameter(1,1.51674e+01);
1362 fFuncWeight->SetParameter(2,4.09941e-01);
96ca9111 1363 fUseWeight=kTRUE;
942a78ce 1364}
1365//_________________________________________________________________________
1366void AliCFTaskVertexingHF::SetPtWeightsFromDataPbPb276overLHC12a17a(){
1367 // ad-hoc weight function from ratio of
1368 // D0 pt spectra in PbPb 2011 0-10% centrality and
1369 // pt spectra from MC production LHC12a17a (PYTHIA Perugia0 with pthard bins)
1370 if(fFuncWeight) delete fFuncWeight;
1371 fFuncWeight=new TF1("funcWeight","[0]+[1]/TMath::Power(x,[2])",0.05,50.);
1372 fFuncWeight->SetParameter(0,1.43116e-02);
1373 fFuncWeight->SetParameter(1,4.37758e+02);
1374 fFuncWeight->SetParameter(2,3.08583);
96ca9111 1375 fUseWeight=kTRUE;
942a78ce 1376}
4fa6d5f8 1377
1d9442de 1378//_________________________________________________________________________
1379void AliCFTaskVertexingHF::SetPtWeightsFromFONLL276overLHC12a17b(){
1380 // weight function from the ratio of the LHC12a17b MC
1381 // and FONLL calculations for pp data
1382 if(fFuncWeight) delete fFuncWeight;
1383 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.);
1384 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);
1385 fUseWeight=kTRUE;
1386}
1387
1388//_________________________________________________________________________
1389void AliCFTaskVertexingHF::SetPtWeightsFromFONLL276andBAMPSoverLHC12a17b(){
1390 // weight function from the ratio of the LHC12a17b MC
1391 // and FONLL calculations for pp data
1392 // corrected by the BAMPS Raa calculation for 30-50% CC
1393 if(fFuncWeight) delete fFuncWeight;
1394 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.);
1395 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);
1396 fUseWeight=kTRUE;
1397}
1398
3825480a 1399//_________________________________________________________________________
1400void AliCFTaskVertexingHF::SetPtWeightsFromFONLL5overLHC13d3(){
1401 // weight function from the ratio of the LHC12a17b MC
1402 // and FONLL calculations for pp data
1403 if(fFuncWeight) delete fFuncWeight;
1404 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.);
1405 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);
1406 fUseWeight=kTRUE;
1407}
1408
ecca15b8 1409//_________________________________________________________________________
1410void AliCFTaskVertexingHF::SetPtWeightsFromFONLL7overLHC10f6a(){
1411 // weight function from the ratio of the LHC12a17b MC
1412 // and FONLL calculations for pp data
1413 if(fFuncWeight) delete fFuncWeight;
1414 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.);
1415 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);
1416 fUseWeight=kTRUE;
1417}
1418
f2703bd2 1419//_________________________________________________________________________
1420Double_t AliCFTaskVertexingHF::GetWeight(Float_t pt)
1421{
5cd139bc 1422 //
1423 // calculating the weight to fill the container
1424 //
f2703bd2 1425
5cd139bc 1426 // FNOLL central:
1427 // p0 = 1.63297e-01 --> 0.322643
1428 // p1 = 2.96275e+00
1429 // p2 = 2.30301e+00
1430 // p3 = 2.50000e+00
1431
1432 // PYTHIA
1433 // p0 = 1.85906e-01 --> 0.36609
1434 // p1 = 1.94635e+00
1435 // p2 = 1.40463e+00
1436 // p3 = 2.50000e+00
1437
1438 Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
1439 Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
1440
1441 Double_t dndpt_func1 = dNdptFit(pt,func1);
1442 if(fUseFlatPtWeight) dndpt_func1 = 1./30.;
1443 Double_t dndpt_func2 = dNdptFit(pt,func2);
1444 AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
1445 return dndpt_func1/dndpt_func2;
f2703bd2 1446}
1447
1448//__________________________________________________________________________________________________
1449Double_t AliCFTaskVertexingHF::dNdptFit(Float_t pt, Double_t* par)
1450{
5cd139bc 1451 //
1452 // calculating dNdpt
1453 //
f2703bd2 1454
5cd139bc 1455 Double_t denom = TMath::Power((pt/par[1]), par[3] );
1456 Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
f2703bd2 1457
5cd139bc 1458 return dNdpt;
f2703bd2 1459}
6c62cb59 1460
17bf1a34 1461//__________________________________________________________________________________________________
1462Double_t AliCFTaskVertexingHF::GetZWeight(Float_t z, Int_t runnumber){
1463 //
1464 // calculating the z-vtx weight for the given run range
1465 //
1466
1467 if(runnumber>146824 || runnumber<146803) return 1.0;
1468
1469 Double_t func1[3] = {1.0, -0.5, 6.5 };
1470 Double_t func2[3] = {1.0, -0.5, 5.5 };
1471
1472 Double_t dzFunc1 = DodzFit(z,func1);
1473 Double_t dzFunc2 = DodzFit(z,func2);
1474
1475 return dzFunc1/dzFunc2;
1476}
1477
1478//__________________________________________________________________________________________________
1479Double_t AliCFTaskVertexingHF::DodzFit(Float_t z, Double_t* par) {
1480
1481 //
1482 // Gaussian z-vtx shape
1483 //
1484 //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
1485
1486 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]);
1487
1488 return value;
1489}
6c62cb59 1490//__________________________________________________________________________________________________
afb24c40 1491Double_t AliCFTaskVertexingHF::GetNchWeight(Int_t nch){
1492 //
1493 // calculates the Nch weight using the measured and generateed Nch distributions
1494 //
1495 if(nch<=0) return 0.;
f26726ed 1496 if(!fHistoMeasNch || !fHistoMCNch) { AliError("Input histos to evaluate Nch weights missing"); return 0.; }
afb24c40 1497 Double_t pMeas=fHistoMeasNch->GetBinContent(fHistoMeasNch->FindBin(nch));
1498 Double_t pMC=fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(nch));
4fa6d5f8 1499 Double_t weight = pMC>0 ? pMeas/pMC : 0.;
1500 if(fUseTrackletsWeight) weight = pMC;
1501 return weight;
afb24c40 1502}
1503//__________________________________________________________________________________________________
1504void AliCFTaskVertexingHF::CreateMeasuredNchHisto(){
1505 // creates historgam with measured multiplcity distribution in pp 7 TeV collisions (from Eur. Phys. J. C (2010) 68: 345–354)
8e1c1209 1506 //
1507 // for Nch > 70 the points were obtained with a double NBD distribution fit
1508 // 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
1509 // fit1->SetParameter(1,1.63); // k parameter
1510 // fit1->SetParameter(2,12.8); // mean multiplicity
1511 //
1512 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 1513 10.50,11.50,12.50,13.50,14.50,15.50,16.50,17.50,18.50,19.50,
1514 20.50,21.50,22.50,23.50,24.50,25.50,26.50,27.50,28.50,29.50,
1515 30.50,31.50,32.50,33.50,34.50,35.50,36.50,37.50,38.50,39.50,
1516 40.50,41.50,42.50,43.50,44.50,45.50,46.50,47.50,48.50,49.50,
1517 50.50,51.50,52.50,53.50,54.50,55.50,56.50,57.50,58.50,59.50,
8e1c1209 1518 60.50,62.50,64.50,66.50,68.50,70.50,72.50,74.50,76.50,78.50,
1519 80.50,82.50,84.50,86.50,88.50,90.50,92.50,94.50,96.50,98.50,
1520 100.50,102.50};
1521 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 1522 0.03359,0.030729,0.028539,0.026575,0.024653,0.0229,0.021325,0.019768,0.018561,0.017187,
1523 0.01604,0.014836,0.013726,0.012576,0.011481,0.010393,0.009502,0.008776,0.008024,0.007452,
1524 0.006851,0.006428,0.00594,0.005515,0.005102,0.00469,0.004162,0.003811,0.003389,0.003071,
1525 0.002708,0.002422,0.002184,0.001968,0.00186,0.00165,0.001577,0.001387,0.001254,0.001118,
1526 0.001037,0.000942,0.000823,0.000736,0.000654,0.000579,0.000512,0.00049,0.00045,0.000355,
8e1c1209 1527 0.000296,0.000265,0.000193,0.00016,0.000126,0.0000851, 0.0000676,0.0000537,0.0000426, 0.0000338,
1528 0.0000268,0.0000213,0.0000166,0.0000133,0.0000106,0.00000837,0.00000662, 0.00000524,0.00000414, 0.00000327,
1529 0.00000258};
afb24c40 1530
1531 if(fHistoMeasNch) delete fHistoMeasNch;
8e1c1209 1532 fHistoMeasNch=new TH1F("hMeaseNch","",81,nchbins);
1533 for(Int_t i=0; i<81; i++){
afb24c40 1534 fHistoMeasNch->SetBinContent(i+1,pch[i]);
1535 fHistoMeasNch->SetBinError(i+1,0.);
1536 }
1537}
1538
5cd139bc 1539//__________________________________________________________________________________________________
6c62cb59 1540Bool_t AliCFTaskVertexingHF::ProcessDs(Int_t recoAnalysisCuts) const{
1541 // processes output of Ds is selected
1542 Bool_t keep=kFALSE;
1543 if(recoAnalysisCuts > 0){
1544 Int_t isKKpi=recoAnalysisCuts&1;
1545 Int_t ispiKK=recoAnalysisCuts&2;
1546 Int_t isPhiKKpi=recoAnalysisCuts&4;
1547 Int_t isPhipiKK=recoAnalysisCuts&8;
1548 Int_t isK0starKKpi=recoAnalysisCuts&16;
1549 Int_t isK0starpiKK=recoAnalysisCuts&32;
1550 if(fDsOption==1){
1551 if(isKKpi && isPhiKKpi) keep=kTRUE;
1552 if(ispiKK && isPhipiKK) keep=kTRUE;
1553 }
1554 else if(fDsOption==2){
1555 if(isKKpi && isK0starKKpi) keep=kTRUE;
1556 if(ispiKK && isK0starpiKK) keep=kTRUE;
1557 }
1558 else if(fDsOption==3)keep=kTRUE;
1559 }
1560 return keep;
1561}
5cd139bc 1562//__________________________________________________________________________________________________
1563Bool_t AliCFTaskVertexingHF::ProcessLctoV0Bachelor(Int_t recoAnalysisCuts) const{
1564 // processes output of Lc->V0+bnachelor is selected
1565 Bool_t keep=kFALSE;
1566 if(recoAnalysisCuts > 0){
1567
1568 Int_t isK0Sp = recoAnalysisCuts&1;
1569 Int_t isLambdaBarpi = recoAnalysisCuts&2;
1570 Int_t isLambdapi = recoAnalysisCuts&4;
1571
1572 if(fLctoV0bachelorOption==1){
1573 if(isK0Sp) keep=kTRUE;
1574 }
1575 else if(fLctoV0bachelorOption==2){
1576 if(isLambdaBarpi) keep=kTRUE;
1577 }
1578 else if(fLctoV0bachelorOption==4){
1579 if(isLambdapi) keep=kTRUE;
1580 }
1581 else if(fLctoV0bachelorOption==7) {
1582 if (isK0Sp || isLambdaBarpi || isLambdapi) keep=kTRUE;
1583 }
1584 }
1585 return keep;
1586}
f26726ed 1587
1588//____________________________________________________________________________
1589TProfile* AliCFTaskVertexingHF::GetEstimatorHistogram(const AliVEvent* event){
1590 // Get Estimator Histogram from period event->GetRunNumber();
1591 //
1592 // If you select SPD tracklets in |eta|<1 you should use type == 1
1593 //
1594
1595 Int_t runNo = event->GetRunNumber();
632f2829 1596 Int_t period = -1; // pp: 0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
1597 // pPb: 0-LHC13b, 1-LHC13c
131192dd 1598
1599 if (fIsPPbData) { // setting run numbers for LHC13 if pPb
1600 if (runNo>195343 && runNo<195484) period = 0;
1601 if (runNo>195528 && runNo<195678) period = 1;
632f2829 1602 if (period<0 || period>1) return 0;
1603 } else { //else assume pp 2010
131192dd 1604 if(runNo>114930 && runNo<117223) period = 0;
1605 if(runNo>119158 && runNo<120830) period = 1;
1606 if(runNo>122373 && runNo<126438) period = 2;
1607 if(runNo>127711 && runNo<130841) period = 3;
1608 if(period<0 || period>3) return 0;
1609 }
632f2829 1610
f26726ed 1611 return fMultEstimatorAvg[period];
1612}