]> 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
Added PbPb default cuts (Rossella)
[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
379592af 378 //loop on the MC event
379
380 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
381 if (!mcArray) {
382 AliError("Could not find Monte-Carlo in AOD");
383 return;
384 }
385 Int_t icountMC = 0;
386 Int_t icountAcc = 0;
387 Int_t icountReco = 0;
388 Int_t icountVertex = 0;
389 Int_t icountRefit = 0;
390 Int_t icountRecoAcc = 0;
391 Int_t icountRecoITSClusters = 0;
392 Int_t icountRecoPPR = 0;
393 Int_t icountRecoPID = 0;
379592af 394 Int_t cquarks = 0;
f2703bd2 395
379592af 396 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
397 if (!mcHeader) {
398 AliError("Could not find MC Header in AOD");
399 return;
400 }
e11ae259 401
402 Double_t* containerInput = new Double_t[fNvar];
403 Double_t* containerInputMC = new Double_t[fNvar];
404
f2703bd2 405
406 AliCFVertexingHF* cfVtxHF=0x0;
407 switch (fDecayChannel){
408 case 2:{
6e2e4f50 409 cfVtxHF = new AliCFVertexingHF2Prong(mcArray, fOriginDselection);
f2703bd2 410 break;
411 }
412 case 21:{
413 // cfVtxHF = new AliCFVertexingHFCascade(mcArray, originDselection); // not there yet
414 break;
415 }
416 case 31:
417 case 32:
418 case 33:{
4e600f58 419 cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel);
f2703bd2 420 break;
421 }
422 case 4:{
423 //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection); // not there yet
424 break;
425 }
426 default:
427 break;
428 }
429
379592af 430 Double_t zPrimVertex = aodVtx ->GetZ();
431 Double_t zMCVertex = mcHeader->GetVtxZ();
f2703bd2 432
379592af 433 //General settings: vertex, feed down and fill reco container with generated values.
434 cfVtxHF->SetRecoPrimVertex(zPrimVertex);
435 cfVtxHF->SetMCPrimaryVertex(zMCVertex);
436 cfVtxHF->SetFillFromGenerated(fFillFromGenerated);
f2703bd2 437
379592af 438 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
379592af 439
f2703bd2 440 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
441
442 // check the MC-level cuts, must be the desidered particle
443 if (!fCFManager->CheckParticleCuts(0, mcPart)) continue; // 0 stands for MC level
444
445 cfVtxHF->SetMCCandidateParam(iPart);
446 cfVtxHF->SetNVar(fNvar);
447 //counting c quarks
448 cquarks += cfVtxHF->MCcquarkCounting(mcPart);
449
6e2e4f50 450 if (!(cfVtxHF->SetLabelArray())){
451 AliDebug(2,Form("Impossible to set the label array (decaychannel = %d)",fDecayChannel));
452 continue;
453 }
454
f2703bd2 455 //check the candiate family at MC level
456 if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) {
6e2e4f50 457 AliDebug(2,Form("Check on the family wrong!!! (decaychannel = %d)",fDecayChannel));
f2703bd2 458 continue;
459 }
460 else{
6e2e4f50 461 AliDebug(2,Form("Check on the family OK!!! (decaychannel = %d)",fDecayChannel));
f2703bd2 462 }
463
464 //Fill the MC container
465 Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
f2703bd2 466 if (mcContainerFilled) {
467 if (fUseWeight)fWeight = GetWeight(containerInputMC[0]);
468 if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) continue;
469 //MC Limited Acceptance
470 if (TMath::Abs(containerInputMC[1]) < 0.5) {
471 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight);
472 AliDebug(3,"MC Lim Acc container filled\n");
473 }
474
475 //MC
476 fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight);
477 icountMC++;
478 AliDebug(3,"MC cointainer filled \n");
479
480 // MC in acceptance
481 // check the MC-Acceptance level cuts
482 // since standard CF functions are not applicable, using Kine Cuts on daughters
483 Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep();
484 if (mcAccepStep){
485 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight);
486 AliDebug(3,"MC acceptance cut passed\n");
487 icountAcc++;
488
489 //MC Vertex step
490 if (fCuts->IsEventSelected(aodEvent)){
491 // filling the container if the vertex is ok
492 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex, fWeight) ;
493 AliDebug(3,"Vertex cut passed and container filled\n");
494 icountVertex++;
495
496 //mc Refit requirement
497 Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, trackCuts);
498 if (mcRefitStep){
499 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight);
500 AliDebug(3,"MC Refit cut passed and container filled\n");
501 icountRefit++;
502 }
503 else{
504 AliDebug(3,"MC Refit cut not passed\n");
505 continue;
506 }
507 }
508 else{
379592af 509 AliDebug (3, "MC vertex step not passed\n");
510 continue;
f2703bd2 511 }
512 }
513 else{
514 AliDebug (3, "MC in acceptance step not passed\n");
515 continue;
516 }
517 }
518 else {
519 AliDebug (3, "MC container not filled\n");
379592af 520 }
521 }
f2703bd2 522
523 if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks));
31da6b05 524 AliDebug(2,Form("Found %i MC particles that are %s!!",icountMC,fPartName.Data()));
525 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Acc cuts!!",icountAcc,fPartName.Data()));
526 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Vertex cuts!!",icountVertex,fPartName.Data()));
527 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Refit cuts!!",icountRefit,fPartName.Data()));
f2703bd2 528
379592af 529 // Now go to rec level
530 fCountMC += icountMC;
531 fCountAcc += icountAcc;
532 fCountVertex+= icountVertex;
533 fCountRefit+= icountRefit;
f2703bd2 534
6e2e4f50 535 AliDebug(2,Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel));
379592af 536
f2703bd2 537 for(Int_t iCandid = 0; iCandid<arrayBranch->GetEntriesFast();iCandid++){
f2703bd2 538 AliAODRecoDecayHF* charmCandidate=0x0;
539 switch (fDecayChannel){
540 case 2:{
541 charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid);
542 break;
543 }
544 case 21:{
545 charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid);
546 break;
547 }
548 case 31:
549 case 32:
550 case 33:{
551 charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid);
552 break;
553 }
554 case 4:{
555 charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid);
556 break;
557 }
558 default:
559 break;
560 }
561
379592af 562 Bool_t unsetvtx=kFALSE;
f2703bd2 563 if(!charmCandidate->GetOwnPrimaryVtx()) {
564 charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
379592af 565 unsetvtx=kTRUE;
566 }
f2703bd2 567
568 Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate);
379592af 569 if (!signAssociation){
f2703bd2 570 charmCandidate = 0x0;
379592af 571 continue;
572 }
3ee5eb83 573
574 Int_t isPartOrAntipart = cfVtxHF->CheckReflexion(fSign);
575 if (isPartOrAntipart == 0){
576 AliDebug(2, Form("The candidate pdg code doesn't match the requirement set in the task (fSign = %d)",fSign));
577 continue;
578 }
579
379592af 580 Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
581 if (recoContFilled){
f2703bd2 582
379592af 583 if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue;
f2703bd2 584
379592af 585 //Reco Step
586 Bool_t recoStep = cfVtxHF->RecoStep();
587 Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
f2703bd2 588
379592af 589 if (recoStep && recoContFilled && vtxCheck){
f2703bd2 590 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
379592af 591 icountReco++;
f2703bd2 592 AliDebug(3,"Reco step passed and container filled\n");
593
379592af 594 //Reco in the acceptance -- take care of UNFOLDING!!!!
595 Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(trackCuts);
596 if (recoAcceptanceStep) {
f2703bd2 597 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
379592af 598 icountRecoAcc++;
f2703bd2 599 AliDebug(3,"Reco acceptance cut passed and container filled\n");
600
379592af 601 if(fAcceptanceUnf){
602 Double_t fill[4]; //fill response matrix
603 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
604 if (bUnfolding) fCorrelation->Fill(fill);
f2703bd2 605 }
606
379592af 607 //Number of ITS cluster requirements
f2703bd2 608 Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
379592af 609 if (recoITSnCluster){
f2703bd2 610 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
379592af 611 icountRecoITSClusters++;
f2703bd2 612 AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
613
614 Bool_t iscutsusingpid = fCuts->GetIsUsePID();
615 Int_t recoAnalysisCuts = -1, recoPidSelection = -1;
616 fCuts->SetUsePID(kFALSE);
617 recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
618
619 if (recoAnalysisCuts > 3){ // Ds case, where more possibilities are considered
620 if (recoAnalysisCuts >= 8){
621 recoAnalysisCuts -= 8; // removing K0star mass
622 }
623 if (recoAnalysisCuts >= 4){
624 recoAnalysisCuts -= 4; // removing Phi mass
625 }
626 }
379592af 627
f2703bd2 628 fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object
629 if (recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart){
630 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight);
631 icountRecoPPR++;
632 AliDebug(3,"Reco Analysis cuts passed and container filled \n");
633 //pid selection
634 //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID);
635 //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){
636 recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
637 if (recoPidSelection == 3 || recoPidSelection == isPartOrAntipart){
638 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight);
639 icountRecoPID++;
640 AliDebug(3,"Reco PID cuts passed and container filled \n");
641 if(!fAcceptanceUnf){
642 Double_t fill[4]; //fill response matrix
643 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
644 if (bUnfolding) fCorrelation->Fill(fill);
645 }
646 }
647 else {
648 AliDebug(3, "Analysis Cuts step not passed \n");
649 continue;
650 }
379592af 651 }
652 else {
f2703bd2 653 AliDebug(3, "PID selection not passed \n");
654 continue;
379592af 655 }
656 }
657 else{
f2703bd2 658 AliDebug(3, "Number of ITS cluster step not passed\n");
659 continue;
379592af 660 }
661 }
662 else{
f2703bd2 663 AliDebug(3, "Reco acceptance step not passed\n");
664 continue;
379592af 665 }
666 }
667 else {
f2703bd2 668 AliDebug(3, "Reco step not passed\n");
669 continue;
379592af 670 }
671 }
672
f2703bd2 673 if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
674 } // end loop on candidate
675
379592af 676 fCountReco+= icountReco;
677 fCountRecoAcc+= icountRecoAcc;
678 fCountRecoITSClusters+= icountRecoITSClusters;
679 fCountRecoPPR+= icountRecoPPR;
680 fCountRecoPID+= icountRecoPID;
681
682 fHistEventsProcessed->Fill(0);
f2703bd2 683
684 delete[] containerInput;
685 containerInput = 0x0;
686 delete[] containerInputMC;
687 containerInputMC = 0x0;
688 delete cfVtxHF;
689
379592af 690}
691
692//___________________________________________________________________________
693void AliCFTaskVertexingHF::Terminate(Option_t*)
694{
695 // The Terminate() function is the last function to be called during
696 // a query. It always runs on the client, it can be used to present
697 // the results graphically or save the results to file.
698
699 AliAnalysisTaskSE::Terminate();
31da6b05 700
701 AliInfo(Form("Found %i MC particles that are %s in MC, in %d events",fCountMC,fPartName.Data(),fEvents));
702 AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, in %d events",fCountAcc,fPartName.Data(),fEvents));
703 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));
704 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));
705 AliInfo(Form("Found %i reco %s that are decaying in %s, in %d events",fCountReco,fPartName.Data(),fDauNames.Data(),fEvents));
706 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));
707 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));
708 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));
709 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 710
711 // draw some example plots....
712
713 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
714 if(!cont) {
715 printf("CONTAINER NOT FOUND\n");
716 return;
717 }
718 // projecting the containers to obtain histograms
719 // first argument = variable, second argument = step
720
31da6b05 721 TH1D* h[3][12];
722 for(Int_t iC=0;iC<12; iC++){
723 // MC-level
724 h[0][iC] = cont->ShowProjection(iC,0);
725 // MC-Acceptance level
726 h[1][iC] = cont->ShowProjection(iC,1);
727 // Reco-level
728 h[2][iC] = cont->ShowProjection(iC,4);
729 }
730 TString titles[12];
731 if(fDecayChannel==31){
732 titles[0]="pT_Dplus (GeV/c)";
733 titles[1]="rapidity";
734 titles[2]="phi (rad)";
735 titles[3]="cT (#mum)";
736 titles[4]="cosPointingAngle";
737 titles[5]="pT_1 (GeV/c)";
738 titles[6]="pT_2 (GeV/c)";
739 titles[7]="pT_3 (GeV/c)";
740 titles[8]="d0_1 (#mum)";
741 titles[9]="d0_2 (#mum)";
742 titles[10]="d0_3 (#mum)";
743 titles[11]="zVertex (cm)";
744 }else{
745 titles[0]="pT_D0 (GeV/c)";
746 titles[1]="rapidity";
747 titles[2]="cosThetaStar";
748 titles[3]="pT_pi (GeV/c)";
749 titles[4]="pT_K (Gev/c)";
750 titles[5]="cT (#mum)";
751 titles[6]="dca (#mum)";
752 titles[7]="d0_pi (#mum)";
753 titles[8]="d0_K (#mum)";
754 titles[9]="d0xd0 (#mum^2)";
755 titles[10]="cosPointingAngle";
756 titles[11]="phi (rad)";
757 }
758 Int_t markers[12]={20,24,21,25,27,28,
759 20,24,21,25,27,28};
760 Int_t colors[3]={2,8,4};
761 for(Int_t iC=0;iC<12; iC++){
762 for(Int_t iStep=0;iStep<3;iStep++){
763 h[iStep][iC]->SetTitle(titles[iC].Data());
764 h[iStep][iC]->GetXaxis()->SetTitle(titles[iC].Data());
765 Double_t maxh=h[iStep][iC]->GetMaximum();
766 h[iStep][iC]->GetYaxis()->SetRangeUser(0,maxh*1.2);
767 h[iStep][iC]->SetMarkerStyle(markers[iC]);
768 h[iStep][iC]->SetMarkerColor(colors[iStep]);
769 }
770 }
771
379592af 772
379592af 773
774 gStyle->SetCanvasColor(0);
775 gStyle->SetFrameFillColor(0);
776 gStyle->SetTitleFillColor(0);
777 gStyle->SetStatColor(0);
778
779 // drawing in 2 separate canvas for a matter of clearity
31da6b05 780 TCanvas * c1 =new TCanvas("c1New","Vars 0,1,2",1100,1600);
379592af 781 c1->Divide(3,3);
31da6b05 782 Int_t iPad=1;
783 for(Int_t iVar=0; iVar<3; iVar++){
784 c1->cd(iPad++);
785 h[0][iVar]->Draw("p");
786 c1->cd(iPad++);
787 h[1][iVar]->Draw("p");
788 c1->cd(iPad++);
789 h[2][iVar]->Draw("p");
790 }
379592af 791
31da6b05 792 TCanvas * c2 =new TCanvas("c2New","Vars 3,4,5",1100,1600);
379592af 793 c2->Divide(3,3);
31da6b05 794 iPad=1;
795 for(Int_t iVar=3; iVar<6; iVar++){
796 c2->cd(iPad++);
797 h[0][iVar]->Draw("p");
798 c2->cd(iPad++);
799 h[1][iVar]->Draw("p");
800 c2->cd(iPad++);
801 h[2][iVar]->Draw("p");
802 }
803
379592af 804
31da6b05 805 TCanvas * c3 =new TCanvas("c3New","Vars 6,7,8",1100,1600);
379592af 806 c3->Divide(3,3);
31da6b05 807 iPad=1;
808 for(Int_t iVar=6; iVar<9; iVar++){
809 c3->cd(iPad++);
810 h[0][iVar]->Draw("p");
811 c3->cd(iPad++);
812 h[1][iVar]->Draw("p");
813 c3->cd(iPad++);
814 h[2][iVar]->Draw("p");
815 }
816
379592af 817
31da6b05 818 TCanvas * c4 =new TCanvas("c4New","Vars 9,10,11",1100,1600);
379592af 819 c4->Divide(3,3);
31da6b05 820 iPad=1;
821 for(Int_t iVar=9; iVar<11; iVar++){
822 c4->cd(iPad++);
823 h[0][iVar]->Draw("p");
824 c4->cd(iPad++);
825 h[1][iVar]->Draw("p");
826 c4->cd(iPad++);
827 h[2][iVar]->Draw("p");
828 }
829
379592af 830
831 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
832
833 TH2D* corr1 =hcorr->Projection(0,2);
834 TH2D* corr2 = hcorr->Projection(1,3);
835
836 TCanvas * c7 =new TCanvas("c7New","",800,400);
837 c7->Divide(2,1);
838 c7->cd(1);
839 corr1->Draw("text");
840 c7->cd(2);
841 corr2->Draw("text");
842
843
844 TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
845
846 corr1->Write();
847 corr2->Write();
31da6b05 848 for(Int_t iC=0;iC<12; iC++){
849 for(Int_t iStep=0;iStep<3;iStep++){
850 h[iStep][iC]->Write(Form("Step%d_%s",iStep,titles[iC].Data()));
851 }
852 }
379592af 853 file_projection->Close();
854
855
856
857 /*
858 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
859 c2->SaveAs("Plots/pTpi_pTK_cT.eps");
860 c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
861 c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
862
863 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
864 c2->SaveAs("Plots/pTpi_pTK_cT.gif");
865 c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
866 c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
867 */
868}
869
870//___________________________________________________________________________
f2703bd2 871void AliCFTaskVertexingHF::UserCreateOutputObjects()
872{
379592af 873 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
874 //TO BE SET BEFORE THE EXECUTION OF THE TASK
875 //
876 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
877
878 //slot #1
879 OpenFile(1);
880 fHistEventsProcessed = new TH1I("CFHFchist0","",1,0,1) ;
881}
882
f2703bd2 883
884//_________________________________________________________________________
885Double_t AliCFTaskVertexingHF::GetWeight(Float_t pt)
886{
887 //
888 // calculating the weight to fill the container
889 //
890
891 // FNOLL central:
892 // p0 = 1.63297e-01 --> 0.322643
893 // p1 = 2.96275e+00
894 // p2 = 2.30301e+00
895 // p3 = 2.50000e+00
896
897 // PYTHIA
898 // p0 = 1.85906e-01 --> 0.36609
899 // p1 = 1.94635e+00
900 // p2 = 1.40463e+00
901 // p3 = 2.50000e+00
902
903 Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
904 Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
905
906 Double_t dndpt_func1 = dNdptFit(pt,func1);
907 Double_t dndpt_func2 = dNdptFit(pt,func2);
908 AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
909 return dndpt_func1/dndpt_func2;
910}
911
912//__________________________________________________________________________________________________
913Double_t AliCFTaskVertexingHF::dNdptFit(Float_t pt, Double_t* par)
914{
915 //
916 // calculating dNdpt
917 //
918
919 Double_t denom = TMath::Power((pt/par[1]), par[3] );
920 Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
921
922 return dNdpt;
923}