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