]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliCFTaskVertexingHF.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bug fix - chack the abs values
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliCFTaskVertexingHF.cxx
CommitLineData
379592af 1/**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//-----------------------------------------------------------------------
17// Class for HF corrections as a function of many variables
18// 6 Steps introduced: MC, MC Acc, Reco, Reco Acc, Reco Acc + ITS Cl,
19// Reco Acc + ITS Cl + PPR cuts
20// 12 variables used: pt, y, cosThetaStar, ptPi, ptK, ct,
21// dca, d0Pi, d0K, d0Pixd0K, cosPointingAngle, phi
22//
23//-----------------------------------------------------------------------
24// Author : C. Zampolli, CERN
25// D. Caffarri, Univ & INFN Padova caffarri@pd.infn.it
26//-----------------------------------------------------------------------
27//-----------------------------------------------------------------------
28// Base class for HF Unfolding (pt and eta)
29// correlation matrix filled at Acceptance and PPR level
30// Author: A.Grelli , Utrecht - agrelli@uu.nl
31//-----------------------------------------------------------------------
32#include <TCanvas.h>
33#include <TParticle.h>
34#include <TDatabasePDG.h>
35#include <TH1I.h>
36#include <TStyle.h>
37#include <TFile.h>
38
39#include "AliCFTaskVertexingHF.h"
40#include "AliStack.h"
41#include "AliMCEvent.h"
42#include "AliCFManager.h"
43#include "AliCFContainer.h"
44#include "AliLog.h"
45#include "AliAnalysisManager.h"
46#include "AliAODHandler.h"
47#include "AliAODEvent.h"
48#include "AliAODRecoDecay.h"
49#include "AliAODRecoDecayHF.h"
50#include "AliAODRecoDecayHF2Prong.h"
f2703bd2 51#include "AliAODRecoDecayHF3Prong.h"
52#include "AliAODRecoDecayHF4Prong.h"
53#include "AliAODRecoCascadeHF.h"
379592af 54#include "AliAODMCParticle.h"
55#include "AliAODMCHeader.h"
56#include "AliESDtrack.h"
57#include "TChain.h"
58#include "THnSparse.h"
59#include "TH2D.h"
60#include "AliESDtrackCuts.h"
f2703bd2 61#include "AliRDHFCuts.h"
379592af 62#include "AliRDHFCutsD0toKpi.h"
f2703bd2 63#include "AliRDHFCutsDplustoKpipi.h"
64#include "AliRDHFCutsDStartoKpipi.h"
65#include "AliRDHFCutsDstoKKpi.h"
66#include "AliRDHFCutsLctopKpi.h"
67#include "AliRDHFCutsD0toKpipipi.h"
379592af 68#include "AliCFVertexingHF2Prong.h"
4e600f58 69#include "AliCFVertexingHF3Prong.h"
367e9aa3 70#include "AliCFVertexingHFCascade.h"
379592af 71#include "AliCFVertexingHF.h"
f2703bd2 72#include "AliAnalysisDataSlot.h"
73#include "AliAnalysisDataContainer.h"
379592af 74
75//__________________________________________________________________________
76AliCFTaskVertexingHF::AliCFTaskVertexingHF() :
f2703bd2 77 AliAnalysisTaskSE(),
f2703bd2 78 fCFManager(0x0),
79 fHistEventsProcessed(0x0),
80 fCorrelation(0x0),
81 fCountMC(0),
82 fCountAcc(0),
83 fCountVertex(0),
84 fCountRefit(0),
85 fCountReco(0),
86 fCountRecoAcc(0),
87 fCountRecoITSClusters(0),
88 fCountRecoPPR(0),
89 fCountRecoPID(0),
90 fEvents(0),
91 fDecayChannel(0),
92 fFillFromGenerated(kFALSE),
93 fOriginDselection(0),
94 fAcceptanceUnf(kTRUE),
95 fCuts(0),
96 fUseWeight(kFALSE),
97 fWeight(1.),
31da6b05 98 fNvar(0),
99 fPartName(""),
3ee5eb83 100 fDauNames(""),
b7af2639 101 fSign(2),
102 fCentralitySelection(kTRUE)
379592af 103{
104 //
105 //Default ctor
106 //
107}
108//___________________________________________________________________________
f2703bd2 109AliCFTaskVertexingHF::AliCFTaskVertexingHF(const Char_t* name, AliRDHFCuts* cuts) :
110 AliAnalysisTaskSE(name),
f2703bd2 111 fCFManager(0x0),
112 fHistEventsProcessed(0x0),
113 fCorrelation(0x0),
114 fCountMC(0),
115 fCountAcc(0),
116 fCountVertex(0),
117 fCountRefit(0),
118 fCountReco(0),
119 fCountRecoAcc(0),
120 fCountRecoITSClusters(0),
121 fCountRecoPPR(0),
122 fCountRecoPID(0),
123 fEvents(0),
124 fDecayChannel(0),
125 fFillFromGenerated(kFALSE),
126 fOriginDselection(0),
127 fAcceptanceUnf(kTRUE),
128 fCuts(cuts),
129 fUseWeight(kFALSE),
130 fWeight(1.),
31da6b05 131 fNvar(0),
132 fPartName(""),
3ee5eb83 133 fDauNames(""),
b7af2639 134 fSign(2),
135 fCentralitySelection(kTRUE)
379592af 136{
137 //
138 // Constructor. Initialization of Inputs and Outputs
139 //
379592af 140 /*
f2703bd2 141 DefineInput(0) and DefineOutput(0)
142 are taken care of by AliAnalysisTaskSE constructor
143 */
379592af 144 DefineOutput(1,TH1I::Class());
145 DefineOutput(2,AliCFContainer::Class());
146 DefineOutput(3,THnSparseD::Class());
f2703bd2 147 DefineOutput(4,AliRDHFCuts::Class());
379592af 148
149 fCuts->PrintAll();
150}
151
152//___________________________________________________________________________
153AliCFTaskVertexingHF& AliCFTaskVertexingHF::operator=(const AliCFTaskVertexingHF& c)
154{
155 //
156 // Assignment operator
157 //
158 if (this!=&c) {
159 AliAnalysisTaskSE::operator=(c) ;
379592af 160 fCFManager = c.fCFManager;
161 fHistEventsProcessed = c.fHistEventsProcessed;
162 fCuts = c.fCuts;
163 }
164 return *this;
165}
166
167//___________________________________________________________________________
168AliCFTaskVertexingHF::AliCFTaskVertexingHF(const AliCFTaskVertexingHF& c) :
f2703bd2 169 AliAnalysisTaskSE(c),
f2703bd2 170 fCFManager(c.fCFManager),
171 fHistEventsProcessed(c.fHistEventsProcessed),
172 fCorrelation(c.fCorrelation),
173 fCountMC(c.fCountMC),
174 fCountAcc(c.fCountAcc),
175 fCountVertex(c.fCountVertex),
176 fCountRefit(c.fCountRefit),
177 fCountReco(c.fCountReco),
178 fCountRecoAcc(c.fCountRecoAcc),
179 fCountRecoITSClusters(c.fCountRecoITSClusters),
180 fCountRecoPPR(c.fCountRecoPPR),
181 fCountRecoPID(c.fCountRecoPID),
182 fEvents(c.fEvents),
183 fDecayChannel(c.fDecayChannel),
184 fFillFromGenerated(c.fFillFromGenerated),
185 fOriginDselection(c.fOriginDselection),
186 fAcceptanceUnf(c.fAcceptanceUnf),
187 fCuts(c.fCuts),
188 fUseWeight(c.fUseWeight),
189 fWeight(c.fWeight),
31da6b05 190 fNvar(c.fNvar),
191 fPartName(c.fPartName),
3ee5eb83 192 fDauNames(c.fDauNames),
b7af2639 193 fSign(c.fSign),
194 fCentralitySelection(c.fCentralitySelection)
379592af 195{
196 //
197 // Copy Constructor
198 //
199}
200
201//___________________________________________________________________________
f2703bd2 202AliCFTaskVertexingHF::~AliCFTaskVertexingHF()
203{
379592af 204 //
205 //destructor
206 //
207 if (fCFManager) delete fCFManager ;
208 if (fHistEventsProcessed) delete fHistEventsProcessed ;
f2703bd2 209 if (fCorrelation) delete fCorrelation ;
379592af 210 if (fCuts) delete fCuts;
211}
212
f2703bd2 213//_________________________________________________________________________-
214void AliCFTaskVertexingHF::Init()
215{
216 //
217 // Initialization
218 //
219
220 if (fDebug>1) printf("AliCFTaskVertexingHF::Init()");
221 AliRDHFCuts *copyfCuts = 0x0;
8a983038 222 if (!fCuts){
223 AliFatal("No cuts defined - Exiting...");
224 return;
225 }
226
f2703bd2 227 switch (fDecayChannel){
228 case 2:{
314ec0fe 229 copyfCuts = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
b7af2639 230 fNvar = 14;
31da6b05 231 fPartName="D0";
232 fDauNames="K+pi";
f2703bd2 233 break;
234 }
235 case 21:{
314ec0fe 236 copyfCuts = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
b7af2639 237 fNvar = 14;
31da6b05 238 fPartName="Dstar";
239 fDauNames="K+pi+pi";
f2703bd2 240 break;
241 }
242 case 31:{
314ec0fe 243 copyfCuts = new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
b7af2639 244 fNvar = 13;
31da6b05 245 fPartName="Dplus";
246 fDauNames="K+pi+pi";
f2703bd2 247 break;
248 }
249 case 32:{
314ec0fe 250 copyfCuts = new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fCuts)));
b7af2639 251 fNvar = 13;
31da6b05 252 fPartName="Lambdac";
253 fDauNames="p+K+pi";
f2703bd2 254 break;
255 }
256 case 33:{
314ec0fe 257 copyfCuts = new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
b7af2639 258 fNvar = 13;
31da6b05 259 fPartName="Ds";
260 fDauNames="K+K+pi";
f2703bd2 261 break;
262 }
263 case 4:{
314ec0fe 264 copyfCuts = new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
b7af2639 265 fNvar = 14;
31da6b05 266 fPartName="D0";
267 fDauNames="K+pi+pi+pi";
f2703bd2 268 break;
269 }
270 default:
271 AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exiting...");
272 break;
273 }
274
275 const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
8a983038 276 if (copyfCuts){
277 copyfCuts->SetName(nameoutput);
278
279 //Post the data
280 PostData(4, copyfCuts);
281 }
282 else{
283 AliFatal("Failing initializing AliRDHFCuts object - Exiting...");
284 }
f2703bd2 285
286 return;
287}
288
379592af 289//_________________________________________________
290void AliCFTaskVertexingHF::UserExec(Option_t *)
291{
379592af 292 //
293 // Main loop function
294 //
295
296 PostData(1,fHistEventsProcessed) ;
297 PostData(2,fCFManager->GetParticleContainer()) ;
298 PostData(3,fCorrelation) ;
299
b7af2639 300
379592af 301 if (fFillFromGenerated){
302 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
303 }
304
305 if (!fInputEvent) {
306 Error("UserExec","NO EVENT FOUND!");
307 return;
308 }
309
310 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
311
f2703bd2 312 TClonesArray *arrayBranch=0;
379592af 313
314 if(!aodEvent && AODEvent() && IsStandardAOD()) {
f2703bd2 315 // In case there is an AOD handler writing a standard AOD, use the AOD
316 // event in memory rather than the input (ESD) event.
317 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
318 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
319 // have to taken from the AOD event hold by the AliAODExtension
320 AliAODHandler* aodHandler = (AliAODHandler*)
321 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
322 if(aodHandler->GetExtensions()) {
323 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
324 AliAODEvent *aodFromExt = ext->GetAOD();
325
326 switch (fDecayChannel){
327 case 2:{
328 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
329 break;
330 }
331 case 21:{
332 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
333 break;
334 }
335 case 31:
336 case 32:
337 case 33:{
338 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
339 break;
340 }
341 case 4:{
342 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
343 break;
344 }
345 default:
346 break;
347 }
348 }
349 }
350 else {
351 switch (fDecayChannel){
352 case 2:{
353 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
354 break;
355 }
356 case 21:{
357 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
358 break;
359 }
360 case 31:
361 case 32:
362 case 33:{
363 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
364 break;
365 }
366 case 4:{
367 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm4Prong");
368 break;
369 }
370 default:
371 break;
372 }
379592af 373 }
374
375 AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
376 if (!aodVtx) return;
f2703bd2 377
378 if (!arrayBranch) {
379 AliError("Could not find array of HF vertices");
380 return;
379592af 381 }
382
383 fEvents++;
6e2e4f50 384
379592af 385 fCFManager->SetRecEventInfo(aodEvent);
386 fCFManager->SetMCEventInfo(aodEvent);
387
379592af 388 //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro.
379592af 389
379592af 390 //loop on the MC event
391
392 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
393 if (!mcArray) {
394 AliError("Could not find Monte-Carlo in AOD");
395 return;
396 }
397 Int_t icountMC = 0;
398 Int_t icountAcc = 0;
399 Int_t icountReco = 0;
400 Int_t icountVertex = 0;
401 Int_t icountRefit = 0;
402 Int_t icountRecoAcc = 0;
403 Int_t icountRecoITSClusters = 0;
404 Int_t icountRecoPPR = 0;
405 Int_t icountRecoPID = 0;
379592af 406 Int_t cquarks = 0;
f2703bd2 407
379592af 408 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
409 if (!mcHeader) {
410 AliError("Could not find MC Header in AOD");
411 return;
412 }
e11ae259 413
414 Double_t* containerInput = new Double_t[fNvar];
415 Double_t* containerInputMC = new Double_t[fNvar];
416
f2703bd2 417
418 AliCFVertexingHF* cfVtxHF=0x0;
419 switch (fDecayChannel){
420 case 2:{
367e9aa3 421 cfVtxHF = new AliCFVertexingHF2Prong(mcArray, fOriginDselection);
422 break;
f2703bd2 423 }
424 case 21:{
367e9aa3 425 cfVtxHF = new AliCFVertexingHFCascade(mcArray, fOriginDselection);
426 break;
f2703bd2 427 }
428 case 31:
429 case 32:
430 case 33:{
4e600f58 431 cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel);
f2703bd2 432 break;
433 }
434 case 4:{
435 //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection); // not there yet
436 break;
437 }
438 default:
439 break;
440 }
8a983038 441 if (!cfVtxHF){
442 AliError("No AliCFVertexingHF initialized");
d49092e3 443 delete[] containerInput;
444 delete[] containerInputMC;
8a983038 445 return;
446 }
f2703bd2 447
379592af 448 Double_t zPrimVertex = aodVtx ->GetZ();
449 Double_t zMCVertex = mcHeader->GetVtxZ();
f2703bd2 450
2bf2e62b 451 AliESDtrackCuts** trackCuts = new AliESDtrackCuts*[cfVtxHF->GetNProngs()];
452 if (fDecayChannel == 21){
453 // for the D*, setting the third element of the array of the track cuts to those for the soft pion
454 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs()-1; iProng++){
455 trackCuts[iProng]=fCuts->GetTrackCuts();
456 }
367e9aa3 457 trackCuts[2] = fCuts->GetTrackCutsSoftPi();
2bf2e62b 458 }
459 else {
460 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs(); iProng++){
461 trackCuts[iProng]=fCuts->GetTrackCuts();
462 }
463 }
464
379592af 465 //General settings: vertex, feed down and fill reco container with generated values.
466 cfVtxHF->SetRecoPrimVertex(zPrimVertex);
467 cfVtxHF->SetMCPrimaryVertex(zMCVertex);
468 cfVtxHF->SetFillFromGenerated(fFillFromGenerated);
b7af2639 469 cfVtxHF->SetNVar(fNvar);
470
471 if (fCentralitySelection)
661d6c29 472 if(fCuts->IsEventSelectedInCentrality(aodEvent)!=0) {
473 delete[] containerInput;
474 delete[] containerInputMC;
475 delete [] trackCuts;
476 return;
477 }
b7af2639 478
479 Float_t centValue = fCuts->GetCentrality(aodEvent);
480 cfVtxHF->SetCentralityValue(centValue);
f2703bd2 481
379592af 482 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
b7af2639 483 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
484 if (!mcPart){
485 AliError("Failed casting particle from MC array!, Skipping particle");
486 continue;
487 }
488 // check the MC-level cuts, must be the desidered particle
489 if (!fCFManager->CheckParticleCuts(0, mcPart)) {
490 continue; // 0 stands for MC level
491 }
492 cfVtxHF->SetMCCandidateParam(iPart);
493
494 //counting c quarks
495 cquarks += cfVtxHF->MCcquarkCounting(mcPart);
496
497 if (!(cfVtxHF->SetLabelArray())){
498 AliDebug(2,Form("Impossible to set the label array (decaychannel = %d)",fDecayChannel));
499 continue;
500 }
2bf2e62b 501
b7af2639 502 //check the candiate family at MC level
503 if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) {
504 AliDebug(2,Form("Check on the family wrong!!! (decaychannel = %d)",fDecayChannel));
505 continue;
506 }
f2703bd2 507 else{
6e2e4f50 508 AliDebug(2,Form("Check on the family OK!!! (decaychannel = %d)",fDecayChannel));
f2703bd2 509 }
510
511 //Fill the MC container
512 Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
f2703bd2 513 if (mcContainerFilled) {
514 if (fUseWeight)fWeight = GetWeight(containerInputMC[0]);
515 if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) continue;
516 //MC Limited Acceptance
517 if (TMath::Abs(containerInputMC[1]) < 0.5) {
518 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight);
519 AliDebug(3,"MC Lim Acc container filled\n");
520 }
521
522 //MC
523 fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight);
524 icountMC++;
525 AliDebug(3,"MC cointainer filled \n");
526
527 // MC in acceptance
528 // check the MC-Acceptance level cuts
529 // since standard CF functions are not applicable, using Kine Cuts on daughters
530 Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep();
531 if (mcAccepStep){
532 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight);
533 AliDebug(3,"MC acceptance cut passed\n");
534 icountAcc++;
535
536 //MC Vertex step
537 if (fCuts->IsEventSelected(aodEvent)){
538 // filling the container if the vertex is ok
539 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex, fWeight) ;
540 AliDebug(3,"Vertex cut passed and container filled\n");
541 icountVertex++;
542
543 //mc Refit requirement
2bf2e62b 544 Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, &trackCuts[0]);
f2703bd2 545 if (mcRefitStep){
546 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight);
547 AliDebug(3,"MC Refit cut passed and container filled\n");
548 icountRefit++;
b7af2639 549
550
f2703bd2 551 }
552 else{
553 AliDebug(3,"MC Refit cut not passed\n");
554 continue;
555 }
556 }
557 else{
379592af 558 AliDebug (3, "MC vertex step not passed\n");
559 continue;
f2703bd2 560 }
561 }
562 else{
563 AliDebug (3, "MC in acceptance step not passed\n");
564 continue;
565 }
566 }
567 else {
568 AliDebug (3, "MC container not filled\n");
379592af 569 }
570 }
f2703bd2 571
572 if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks));
31da6b05 573 AliDebug(2,Form("Found %i MC particles that are %s!!",icountMC,fPartName.Data()));
574 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Acc cuts!!",icountAcc,fPartName.Data()));
575 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Vertex cuts!!",icountVertex,fPartName.Data()));
576 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Refit cuts!!",icountRefit,fPartName.Data()));
f2703bd2 577
379592af 578 // Now go to rec level
579 fCountMC += icountMC;
580 fCountAcc += icountAcc;
581 fCountVertex+= icountVertex;
582 fCountRefit+= icountRefit;
f2703bd2 583
6e2e4f50 584 AliDebug(2,Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel));
379592af 585
f2703bd2 586 for(Int_t iCandid = 0; iCandid<arrayBranch->GetEntriesFast();iCandid++){
f2703bd2 587 AliAODRecoDecayHF* charmCandidate=0x0;
588 switch (fDecayChannel){
589 case 2:{
590 charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid);
591 break;
592 }
593 case 21:{
594 charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid);
595 break;
596 }
597 case 31:
598 case 32:
599 case 33:{
600 charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid);
601 break;
602 }
603 case 4:{
604 charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid);
605 break;
606 }
607 default:
608 break;
609 }
610
379592af 611 Bool_t unsetvtx=kFALSE;
f2703bd2 612 if(!charmCandidate->GetOwnPrimaryVtx()) {
613 charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
379592af 614 unsetvtx=kTRUE;
615 }
f2703bd2 616
617 Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate);
379592af 618 if (!signAssociation){
f2703bd2 619 charmCandidate = 0x0;
379592af 620 continue;
621 }
3ee5eb83 622
623 Int_t isPartOrAntipart = cfVtxHF->CheckReflexion(fSign);
624 if (isPartOrAntipart == 0){
625 AliDebug(2, Form("The candidate pdg code doesn't match the requirement set in the task (fSign = %d)",fSign));
626 continue;
627 }
628
379592af 629 Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
630 if (recoContFilled){
f2703bd2 631
379592af 632 if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue;
f2703bd2 633
379592af 634 //Reco Step
635 Bool_t recoStep = cfVtxHF->RecoStep();
636 Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
f2703bd2 637
379592af 638 if (recoStep && recoContFilled && vtxCheck){
f2703bd2 639 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
379592af 640 icountReco++;
f2703bd2 641 AliDebug(3,"Reco step passed and container filled\n");
642
379592af 643 //Reco in the acceptance -- take care of UNFOLDING!!!!
2bf2e62b 644 Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(&trackCuts[0]);
379592af 645 if (recoAcceptanceStep) {
f2703bd2 646 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
379592af 647 icountRecoAcc++;
f2703bd2 648 AliDebug(3,"Reco acceptance cut passed and container filled\n");
649
379592af 650 if(fAcceptanceUnf){
651 Double_t fill[4]; //fill response matrix
652 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
653 if (bUnfolding) fCorrelation->Fill(fill);
f2703bd2 654 }
655
379592af 656 //Number of ITS cluster requirements
f2703bd2 657 Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
379592af 658 if (recoITSnCluster){
f2703bd2 659 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
379592af 660 icountRecoITSClusters++;
f2703bd2 661 AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
662
663 Bool_t iscutsusingpid = fCuts->GetIsUsePID();
664 Int_t recoAnalysisCuts = -1, recoPidSelection = -1;
665 fCuts->SetUsePID(kFALSE);
666 recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
667
668 if (recoAnalysisCuts > 3){ // Ds case, where more possibilities are considered
669 if (recoAnalysisCuts >= 8){
670 recoAnalysisCuts -= 8; // removing K0star mass
671 }
672 if (recoAnalysisCuts >= 4){
673 recoAnalysisCuts -= 4; // removing Phi mass
674 }
675 }
379592af 676
f2703bd2 677 fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object
678 if (recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart){
679 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight);
680 icountRecoPPR++;
681 AliDebug(3,"Reco Analysis cuts passed and container filled \n");
682 //pid selection
683 //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID);
684 //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){
685 recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
686 if (recoPidSelection == 3 || recoPidSelection == isPartOrAntipart){
687 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight);
688 icountRecoPID++;
689 AliDebug(3,"Reco PID cuts passed and container filled \n");
690 if(!fAcceptanceUnf){
691 Double_t fill[4]; //fill response matrix
692 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
693 if (bUnfolding) fCorrelation->Fill(fill);
694 }
695 }
696 else {
697 AliDebug(3, "Analysis Cuts step not passed \n");
698 continue;
699 }
379592af 700 }
701 else {
f2703bd2 702 AliDebug(3, "PID selection not passed \n");
703 continue;
379592af 704 }
705 }
706 else{
f2703bd2 707 AliDebug(3, "Number of ITS cluster step not passed\n");
708 continue;
379592af 709 }
710 }
711 else{
f2703bd2 712 AliDebug(3, "Reco acceptance step not passed\n");
713 continue;
379592af 714 }
715 }
716 else {
f2703bd2 717 AliDebug(3, "Reco step not passed\n");
718 continue;
379592af 719 }
720 }
721
f2703bd2 722 if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
723 } // end loop on candidate
724
379592af 725 fCountReco+= icountReco;
726 fCountRecoAcc+= icountRecoAcc;
727 fCountRecoITSClusters+= icountRecoITSClusters;
728 fCountRecoPPR+= icountRecoPPR;
729 fCountRecoPID+= icountRecoPID;
730
731 fHistEventsProcessed->Fill(0);
f2703bd2 732
733 delete[] containerInput;
f2703bd2 734 delete[] containerInputMC;
f2703bd2 735 delete cfVtxHF;
2bf2e62b 736 if (trackCuts){
737 // for (Int_t i=0; i<cfVtxHF->GetNProngs(); i++){
738 // delete [] trackCuts[i];
739 // }
740 delete [] trackCuts;
741 }
379592af 742}
743
744//___________________________________________________________________________
745void AliCFTaskVertexingHF::Terminate(Option_t*)
746{
747 // The Terminate() function is the last function to be called during
748 // a query. It always runs on the client, it can be used to present
749 // the results graphically or save the results to file.
750
751 AliAnalysisTaskSE::Terminate();
31da6b05 752
753 AliInfo(Form("Found %i MC particles that are %s in MC, in %d events",fCountMC,fPartName.Data(),fEvents));
754 AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, in %d events",fCountAcc,fPartName.Data(),fEvents));
755 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));
756 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));
757 AliInfo(Form("Found %i reco %s that are decaying in %s, in %d events",fCountReco,fPartName.Data(),fDauNames.Data(),fEvents));
758 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));
759 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));
760 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));
761 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 762
763 // draw some example plots....
764
765 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
766 if(!cont) {
767 printf("CONTAINER NOT FOUND\n");
768 return;
769 }
770 // projecting the containers to obtain histograms
771 // first argument = variable, second argument = step
772
31da6b05 773 TH1D* h[3][12];
774 for(Int_t iC=0;iC<12; iC++){
775 // MC-level
776 h[0][iC] = cont->ShowProjection(iC,0);
777 // MC-Acceptance level
778 h[1][iC] = cont->ShowProjection(iC,1);
779 // Reco-level
780 h[2][iC] = cont->ShowProjection(iC,4);
781 }
782 TString titles[12];
783 if(fDecayChannel==31){
784 titles[0]="pT_Dplus (GeV/c)";
785 titles[1]="rapidity";
786 titles[2]="phi (rad)";
787 titles[3]="cT (#mum)";
788 titles[4]="cosPointingAngle";
789 titles[5]="pT_1 (GeV/c)";
790 titles[6]="pT_2 (GeV/c)";
791 titles[7]="pT_3 (GeV/c)";
792 titles[8]="d0_1 (#mum)";
793 titles[9]="d0_2 (#mum)";
794 titles[10]="d0_3 (#mum)";
795 titles[11]="zVertex (cm)";
796 }else{
797 titles[0]="pT_D0 (GeV/c)";
798 titles[1]="rapidity";
799 titles[2]="cosThetaStar";
800 titles[3]="pT_pi (GeV/c)";
801 titles[4]="pT_K (Gev/c)";
802 titles[5]="cT (#mum)";
803 titles[6]="dca (#mum)";
804 titles[7]="d0_pi (#mum)";
805 titles[8]="d0_K (#mum)";
806 titles[9]="d0xd0 (#mum^2)";
807 titles[10]="cosPointingAngle";
808 titles[11]="phi (rad)";
809 }
810 Int_t markers[12]={20,24,21,25,27,28,
811 20,24,21,25,27,28};
812 Int_t colors[3]={2,8,4};
813 for(Int_t iC=0;iC<12; iC++){
814 for(Int_t iStep=0;iStep<3;iStep++){
815 h[iStep][iC]->SetTitle(titles[iC].Data());
816 h[iStep][iC]->GetXaxis()->SetTitle(titles[iC].Data());
817 Double_t maxh=h[iStep][iC]->GetMaximum();
818 h[iStep][iC]->GetYaxis()->SetRangeUser(0,maxh*1.2);
819 h[iStep][iC]->SetMarkerStyle(markers[iC]);
820 h[iStep][iC]->SetMarkerColor(colors[iStep]);
821 }
822 }
823
379592af 824
379592af 825
826 gStyle->SetCanvasColor(0);
827 gStyle->SetFrameFillColor(0);
828 gStyle->SetTitleFillColor(0);
829 gStyle->SetStatColor(0);
830
831 // drawing in 2 separate canvas for a matter of clearity
31da6b05 832 TCanvas * c1 =new TCanvas("c1New","Vars 0,1,2",1100,1600);
379592af 833 c1->Divide(3,3);
31da6b05 834 Int_t iPad=1;
835 for(Int_t iVar=0; iVar<3; iVar++){
836 c1->cd(iPad++);
837 h[0][iVar]->Draw("p");
838 c1->cd(iPad++);
839 h[1][iVar]->Draw("p");
840 c1->cd(iPad++);
841 h[2][iVar]->Draw("p");
842 }
379592af 843
31da6b05 844 TCanvas * c2 =new TCanvas("c2New","Vars 3,4,5",1100,1600);
379592af 845 c2->Divide(3,3);
31da6b05 846 iPad=1;
847 for(Int_t iVar=3; iVar<6; iVar++){
848 c2->cd(iPad++);
849 h[0][iVar]->Draw("p");
850 c2->cd(iPad++);
851 h[1][iVar]->Draw("p");
852 c2->cd(iPad++);
853 h[2][iVar]->Draw("p");
854 }
855
379592af 856
31da6b05 857 TCanvas * c3 =new TCanvas("c3New","Vars 6,7,8",1100,1600);
379592af 858 c3->Divide(3,3);
31da6b05 859 iPad=1;
860 for(Int_t iVar=6; iVar<9; iVar++){
861 c3->cd(iPad++);
862 h[0][iVar]->Draw("p");
863 c3->cd(iPad++);
864 h[1][iVar]->Draw("p");
865 c3->cd(iPad++);
866 h[2][iVar]->Draw("p");
867 }
868
379592af 869
31da6b05 870 TCanvas * c4 =new TCanvas("c4New","Vars 9,10,11",1100,1600);
379592af 871 c4->Divide(3,3);
31da6b05 872 iPad=1;
873 for(Int_t iVar=9; iVar<11; iVar++){
874 c4->cd(iPad++);
875 h[0][iVar]->Draw("p");
876 c4->cd(iPad++);
877 h[1][iVar]->Draw("p");
878 c4->cd(iPad++);
879 h[2][iVar]->Draw("p");
880 }
881
379592af 882
883 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
884
885 TH2D* corr1 =hcorr->Projection(0,2);
886 TH2D* corr2 = hcorr->Projection(1,3);
887
888 TCanvas * c7 =new TCanvas("c7New","",800,400);
889 c7->Divide(2,1);
890 c7->cd(1);
891 corr1->Draw("text");
892 c7->cd(2);
893 corr2->Draw("text");
894
895
896 TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
897
898 corr1->Write();
899 corr2->Write();
31da6b05 900 for(Int_t iC=0;iC<12; iC++){
901 for(Int_t iStep=0;iStep<3;iStep++){
902 h[iStep][iC]->Write(Form("Step%d_%s",iStep,titles[iC].Data()));
903 }
904 }
379592af 905 file_projection->Close();
906
907
908
909 /*
910 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
911 c2->SaveAs("Plots/pTpi_pTK_cT.eps");
912 c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
913 c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
914
915 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
916 c2->SaveAs("Plots/pTpi_pTK_cT.gif");
917 c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
918 c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
919 */
920}
921
922//___________________________________________________________________________
f2703bd2 923void AliCFTaskVertexingHF::UserCreateOutputObjects()
924{
379592af 925 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
926 //TO BE SET BEFORE THE EXECUTION OF THE TASK
927 //
928 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
929
930 //slot #1
931 OpenFile(1);
932 fHistEventsProcessed = new TH1I("CFHFchist0","",1,0,1) ;
933}
934
f2703bd2 935
936//_________________________________________________________________________
937Double_t AliCFTaskVertexingHF::GetWeight(Float_t pt)
938{
939 //
940 // calculating the weight to fill the container
941 //
942
943 // FNOLL central:
944 // p0 = 1.63297e-01 --> 0.322643
945 // p1 = 2.96275e+00
946 // p2 = 2.30301e+00
947 // p3 = 2.50000e+00
948
949 // PYTHIA
950 // p0 = 1.85906e-01 --> 0.36609
951 // p1 = 1.94635e+00
952 // p2 = 1.40463e+00
953 // p3 = 2.50000e+00
954
955 Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
956 Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
957
958 Double_t dndpt_func1 = dNdptFit(pt,func1);
959 Double_t dndpt_func2 = dNdptFit(pt,func2);
960 AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
961 return dndpt_func1/dndpt_func2;
962}
963
964//__________________________________________________________________________________________________
965Double_t AliCFTaskVertexingHF::dNdptFit(Float_t pt, Double_t* par)
966{
967 //
968 // calculating dNdpt
969 //
970
971 Double_t denom = TMath::Power((pt/par[1]), par[3] );
972 Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
973
974 return dNdpt;
975}