]> 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
Move AliError and AliWarning prints to AliDebug
[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"
f2703bd2 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(),
77 fPDG(0),
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.),
98 fNvar(0)
379592af 99{
100 //
101 //Default ctor
102 //
103}
104//___________________________________________________________________________
f2703bd2 105AliCFTaskVertexingHF::AliCFTaskVertexingHF(const Char_t* name, AliRDHFCuts* cuts) :
106 AliAnalysisTaskSE(name),
107 fPDG(0),
108 fCFManager(0x0),
109 fHistEventsProcessed(0x0),
110 fCorrelation(0x0),
111 fCountMC(0),
112 fCountAcc(0),
113 fCountVertex(0),
114 fCountRefit(0),
115 fCountReco(0),
116 fCountRecoAcc(0),
117 fCountRecoITSClusters(0),
118 fCountRecoPPR(0),
119 fCountRecoPID(0),
120 fEvents(0),
121 fDecayChannel(0),
122 fFillFromGenerated(kFALSE),
123 fOriginDselection(0),
124 fAcceptanceUnf(kTRUE),
125 fCuts(cuts),
126 fUseWeight(kFALSE),
127 fWeight(1.),
128 fNvar(0)
379592af 129{
130 //
131 // Constructor. Initialization of Inputs and Outputs
132 //
379592af 133 /*
f2703bd2 134 DefineInput(0) and DefineOutput(0)
135 are taken care of by AliAnalysisTaskSE constructor
136 */
379592af 137 DefineOutput(1,TH1I::Class());
138 DefineOutput(2,AliCFContainer::Class());
139 DefineOutput(3,THnSparseD::Class());
f2703bd2 140 DefineOutput(4,AliRDHFCuts::Class());
379592af 141
142 fCuts->PrintAll();
143}
144
145//___________________________________________________________________________
146AliCFTaskVertexingHF& AliCFTaskVertexingHF::operator=(const AliCFTaskVertexingHF& c)
147{
148 //
149 // Assignment operator
150 //
151 if (this!=&c) {
152 AliAnalysisTaskSE::operator=(c) ;
379592af 153 fCFManager = c.fCFManager;
154 fHistEventsProcessed = c.fHistEventsProcessed;
155 fCuts = c.fCuts;
156 }
157 return *this;
158}
159
160//___________________________________________________________________________
161AliCFTaskVertexingHF::AliCFTaskVertexingHF(const AliCFTaskVertexingHF& c) :
f2703bd2 162 AliAnalysisTaskSE(c),
163 fPDG(c.fPDG),
164 fCFManager(c.fCFManager),
165 fHistEventsProcessed(c.fHistEventsProcessed),
166 fCorrelation(c.fCorrelation),
167 fCountMC(c.fCountMC),
168 fCountAcc(c.fCountAcc),
169 fCountVertex(c.fCountVertex),
170 fCountRefit(c.fCountRefit),
171 fCountReco(c.fCountReco),
172 fCountRecoAcc(c.fCountRecoAcc),
173 fCountRecoITSClusters(c.fCountRecoITSClusters),
174 fCountRecoPPR(c.fCountRecoPPR),
175 fCountRecoPID(c.fCountRecoPID),
176 fEvents(c.fEvents),
177 fDecayChannel(c.fDecayChannel),
178 fFillFromGenerated(c.fFillFromGenerated),
179 fOriginDselection(c.fOriginDselection),
180 fAcceptanceUnf(c.fAcceptanceUnf),
181 fCuts(c.fCuts),
182 fUseWeight(c.fUseWeight),
183 fWeight(c.fWeight),
184 fNvar(c.fNvar)
379592af 185{
186 //
187 // Copy Constructor
188 //
189}
190
191//___________________________________________________________________________
f2703bd2 192AliCFTaskVertexingHF::~AliCFTaskVertexingHF()
193{
379592af 194 //
195 //destructor
196 //
197 if (fCFManager) delete fCFManager ;
198 if (fHistEventsProcessed) delete fHistEventsProcessed ;
f2703bd2 199 if (fCorrelation) delete fCorrelation ;
379592af 200 if (fCuts) delete fCuts;
201}
202
f2703bd2 203//_________________________________________________________________________-
204void AliCFTaskVertexingHF::Init()
205{
206 //
207 // Initialization
208 //
209
210 if (fDebug>1) printf("AliCFTaskVertexingHF::Init()");
211 AliRDHFCuts *copyfCuts = 0x0;
212
213
214 switch (fDecayChannel){
215 case 2:{
216 copyfCuts = new AliRDHFCutsD0toKpi(*(dynamic_cast<AliRDHFCutsD0toKpi*>(fCuts)));
217 fNvar = 13;
218 break;
219 }
220 case 21:{
221 copyfCuts = new AliRDHFCutsDStartoKpipi(*(dynamic_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
222 fNvar = 13;
223 break;
224 }
225 case 31:{
226 copyfCuts = new AliRDHFCutsDplustoKpipi(*(dynamic_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
227 fNvar = 12;
228 break;
229 }
230 case 32:{
231 copyfCuts = new AliRDHFCutsLctopKpi(*(dynamic_cast<AliRDHFCutsLctopKpi*>(fCuts)));
232 fNvar = 13;
233 break;
234 }
235 case 33:{
236 copyfCuts = new AliRDHFCutsDstoKKpi(*(dynamic_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
237 fNvar = 13;
238 break;
239 }
240 case 4:{
241 copyfCuts = new AliRDHFCutsD0toKpipipi(*(dynamic_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
242 fNvar = 13;
243 break;
244 }
245 default:
246 AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exiting...");
247 break;
248 }
249
250 const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
251 copyfCuts->SetName(nameoutput);
252
253 //Post the data
254 PostData(4, copyfCuts);
255
256 return;
257}
258
379592af 259//_________________________________________________
260void AliCFTaskVertexingHF::UserExec(Option_t *)
261{
379592af 262 //
263 // Main loop function
264 //
265
266 PostData(1,fHistEventsProcessed) ;
267 PostData(2,fCFManager->GetParticleContainer()) ;
268 PostData(3,fCorrelation) ;
269
270 AliESDtrackCuts* trackCuts = fCuts->GetTrackCuts();
271
272 if (fFillFromGenerated){
273 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
274 }
275
276 if (!fInputEvent) {
277 Error("UserExec","NO EVENT FOUND!");
278 return;
279 }
280
281 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
282
f2703bd2 283 TClonesArray *arrayBranch=0;
379592af 284
285 if(!aodEvent && AODEvent() && IsStandardAOD()) {
f2703bd2 286 // In case there is an AOD handler writing a standard AOD, use the AOD
287 // event in memory rather than the input (ESD) event.
288 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
289 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
290 // have to taken from the AOD event hold by the AliAODExtension
291 AliAODHandler* aodHandler = (AliAODHandler*)
292 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
293 if(aodHandler->GetExtensions()) {
294 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
295 AliAODEvent *aodFromExt = ext->GetAOD();
296
297 switch (fDecayChannel){
298 case 2:{
299 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
300 break;
301 }
302 case 21:{
303 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
304 break;
305 }
306 case 31:
307 case 32:
308 case 33:{
309 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
310 break;
311 }
312 case 4:{
313 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
314 break;
315 }
316 default:
317 break;
318 }
319 }
320 }
321 else {
322 switch (fDecayChannel){
323 case 2:{
324 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
325 break;
326 }
327 case 21:{
328 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
329 break;
330 }
331 case 31:
332 case 32:
333 case 33:{
334 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
335 break;
336 }
337 case 4:{
338 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm4Prong");
339 break;
340 }
341 default:
342 break;
343 }
379592af 344 }
345
346 AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
347 if (!aodVtx) return;
f2703bd2 348
349 if (!arrayBranch) {
350 AliError("Could not find array of HF vertices");
351 return;
379592af 352 }
353
354 fEvents++;
355 if (fEvents%10000 == 0) AliDebug(2,Form("Event %d",fEvents));
356
357 fCFManager->SetRecEventInfo(aodEvent);
358 fCFManager->SetMCEventInfo(aodEvent);
359
379592af 360 //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro.
f2703bd2 361 // Int_t nVar = 13;
379592af 362
f2703bd2 363 Double_t* containerInput = new Double_t[fNvar];
364 Double_t* containerInputMC = new Double_t[fNvar];
379592af 365
366 //loop on the MC event
367
368 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
369 if (!mcArray) {
370 AliError("Could not find Monte-Carlo in AOD");
371 return;
372 }
373 Int_t icountMC = 0;
374 Int_t icountAcc = 0;
375 Int_t icountReco = 0;
376 Int_t icountVertex = 0;
377 Int_t icountRefit = 0;
378 Int_t icountRecoAcc = 0;
379 Int_t icountRecoITSClusters = 0;
380 Int_t icountRecoPPR = 0;
381 Int_t icountRecoPID = 0;
379592af 382 Int_t cquarks = 0;
379592af 383 UShort_t originDselection = 0;
f2703bd2 384
379592af 385 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
386 if (!mcHeader) {
387 AliError("Could not find MC Header in AOD");
388 return;
389 }
f2703bd2 390
391 AliCFVertexingHF* cfVtxHF=0x0;
392 switch (fDecayChannel){
393 case 2:{
394 cfVtxHF = new AliCFVertexingHF2Prong(mcArray, originDselection);
395 break;
396 }
397 case 21:{
398 // cfVtxHF = new AliCFVertexingHFCascade(mcArray, originDselection); // not there yet
399 break;
400 }
401 case 31:
402 case 32:
403 case 33:{
404 //cfVtxHF = new AliCFVertexingHF3Prong(mcArray, originDselection, fDecayChannel);
405 break;
406 }
407 case 4:{
408 //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection); // not there yet
409 break;
410 }
411 default:
412 break;
413 }
414
379592af 415 Double_t zPrimVertex = aodVtx ->GetZ();
416 Double_t zMCVertex = mcHeader->GetVtxZ();
f2703bd2 417
379592af 418 //General settings: vertex, feed down and fill reco container with generated values.
419 cfVtxHF->SetRecoPrimVertex(zPrimVertex);
420 cfVtxHF->SetMCPrimaryVertex(zMCVertex);
421 cfVtxHF->SetFillFromGenerated(fFillFromGenerated);
f2703bd2 422
379592af 423 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
379592af 424
f2703bd2 425 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
426
427 // check the MC-level cuts, must be the desidered particle
428 if (!fCFManager->CheckParticleCuts(0, mcPart)) continue; // 0 stands for MC level
429
430 cfVtxHF->SetMCCandidateParam(iPart);
431 cfVtxHF->SetNVar(fNvar);
432 //counting c quarks
433 cquarks += cfVtxHF->MCcquarkCounting(mcPart);
434
435 //check the candiate family at MC level
436 if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) {
437 Printf("Check on the family wrong!!! (decaychannel = %d)",fDecayChannel);
438 continue;
439 }
440 else{
441 Printf("Check on the family OK!!! (decaychannel = %d)",fDecayChannel);
442 }
443
444 //Fill the MC container
445 Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
446
447 if (mcContainerFilled) {
448 if (fUseWeight)fWeight = GetWeight(containerInputMC[0]);
449 if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) continue;
450 //MC Limited Acceptance
451 if (TMath::Abs(containerInputMC[1]) < 0.5) {
452 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight);
453 AliDebug(3,"MC Lim Acc container filled\n");
454 }
455
456 //MC
457 fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight);
458 icountMC++;
459 AliDebug(3,"MC cointainer filled \n");
460
461 // MC in acceptance
462 // check the MC-Acceptance level cuts
463 // since standard CF functions are not applicable, using Kine Cuts on daughters
464 Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep();
465 if (mcAccepStep){
466 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight);
467 AliDebug(3,"MC acceptance cut passed\n");
468 icountAcc++;
469
470 //MC Vertex step
471 if (fCuts->IsEventSelected(aodEvent)){
472 // filling the container if the vertex is ok
473 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex, fWeight) ;
474 AliDebug(3,"Vertex cut passed and container filled\n");
475 icountVertex++;
476
477 //mc Refit requirement
478 Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, trackCuts);
479 if (mcRefitStep){
480 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight);
481 AliDebug(3,"MC Refit cut passed and container filled\n");
482 icountRefit++;
483 }
484 else{
485 AliDebug(3,"MC Refit cut not passed\n");
486 continue;
487 }
488 }
489 else{
379592af 490 AliDebug (3, "MC vertex step not passed\n");
491 continue;
f2703bd2 492 }
493 }
494 else{
495 AliDebug (3, "MC in acceptance step not passed\n");
496 continue;
497 }
498 }
499 else {
500 AliDebug (3, "MC container not filled\n");
379592af 501 }
502 }
f2703bd2 503
504 if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks));
379592af 505 AliDebug(2,Form("Found %i MC particles that are D0!!",icountMC));
506 AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Acc cuts!!",icountAcc));
507 AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Vertex cuts!!",icountVertex));
508 AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Refit cuts!!",icountRefit));
f2703bd2 509
379592af 510 // Now go to rec level
511 fCountMC += icountMC;
512 fCountAcc += icountAcc;
513 fCountVertex+= icountVertex;
514 fCountRefit+= icountRefit;
f2703bd2 515
516 AliDebug(2,Form("Found %d vertices",arrayBranch->GetEntriesFast()));
517 AliInfo(Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel));
379592af 518
f2703bd2 519 for(Int_t iCandid = 0; iCandid<arrayBranch->GetEntriesFast();iCandid++){
520 Printf("iCandid = %d", iCandid);
379592af 521
f2703bd2 522 AliAODRecoDecayHF* charmCandidate=0x0;
523 switch (fDecayChannel){
524 case 2:{
525 charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid);
526 break;
527 }
528 case 21:{
529 charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid);
530 break;
531 }
532 case 31:
533 case 32:
534 case 33:{
535 charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid);
536 break;
537 }
538 case 4:{
539 charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid);
540 break;
541 }
542 default:
543 break;
544 }
545
379592af 546 Bool_t unsetvtx=kFALSE;
f2703bd2 547 if(!charmCandidate->GetOwnPrimaryVtx()) {
548 charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
379592af 549 unsetvtx=kTRUE;
550 }
f2703bd2 551
552 Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate);
379592af 553 if (!signAssociation){
f2703bd2 554 charmCandidate = 0x0;
379592af 555 continue;
556 }
f2703bd2 557
558 Int_t isPartOrAntipart = cfVtxHF->CheckReflexion();
559
379592af 560 Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
561 if (recoContFilled){
f2703bd2 562
379592af 563 if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue;
f2703bd2 564
379592af 565 //Reco Step
566 Bool_t recoStep = cfVtxHF->RecoStep();
567 Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
f2703bd2 568
379592af 569 if (recoStep && recoContFilled && vtxCheck){
f2703bd2 570 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
379592af 571 icountReco++;
f2703bd2 572 AliDebug(3,"Reco step passed and container filled\n");
573
379592af 574 //Reco in the acceptance -- take care of UNFOLDING!!!!
575 Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(trackCuts);
576 if (recoAcceptanceStep) {
f2703bd2 577 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
379592af 578 icountRecoAcc++;
f2703bd2 579 AliDebug(3,"Reco acceptance cut passed and container filled\n");
580
379592af 581 if(fAcceptanceUnf){
582 Double_t fill[4]; //fill response matrix
583 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
584 if (bUnfolding) fCorrelation->Fill(fill);
f2703bd2 585 }
586
379592af 587 //Number of ITS cluster requirements
f2703bd2 588 Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
379592af 589 if (recoITSnCluster){
f2703bd2 590 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
379592af 591 icountRecoITSClusters++;
f2703bd2 592 AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
593
594 Bool_t iscutsusingpid = fCuts->GetIsUsePID();
595 Int_t recoAnalysisCuts = -1, recoPidSelection = -1;
596 fCuts->SetUsePID(kFALSE);
597 recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
598
599 if (recoAnalysisCuts > 3){ // Ds case, where more possibilities are considered
600 if (recoAnalysisCuts >= 8){
601 recoAnalysisCuts -= 8; // removing K0star mass
602 }
603 if (recoAnalysisCuts >= 4){
604 recoAnalysisCuts -= 4; // removing Phi mass
605 }
606 }
379592af 607
f2703bd2 608 fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object
609 if (recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart){
610 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight);
611 icountRecoPPR++;
612 AliDebug(3,"Reco Analysis cuts passed and container filled \n");
613 //pid selection
614 //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID);
615 //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){
616 recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
617 if (recoPidSelection == 3 || recoPidSelection == isPartOrAntipart){
618 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight);
619 icountRecoPID++;
620 AliDebug(3,"Reco PID cuts passed and container filled \n");
621 if(!fAcceptanceUnf){
622 Double_t fill[4]; //fill response matrix
623 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
624 if (bUnfolding) fCorrelation->Fill(fill);
625 }
626 }
627 else {
628 AliDebug(3, "Analysis Cuts step not passed \n");
629 continue;
630 }
379592af 631 }
632 else {
f2703bd2 633 AliDebug(3, "PID selection not passed \n");
634 continue;
379592af 635 }
636 }
637 else{
f2703bd2 638 AliDebug(3, "Number of ITS cluster step not passed\n");
639 continue;
379592af 640 }
641 }
642 else{
f2703bd2 643 AliDebug(3, "Reco acceptance step not passed\n");
644 continue;
379592af 645 }
646 }
647 else {
f2703bd2 648 AliDebug(3, "Reco step not passed\n");
649 continue;
379592af 650 }
651 }
652
f2703bd2 653 if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
654 } // end loop on candidate
655
379592af 656 fCountReco+= icountReco;
657 fCountRecoAcc+= icountRecoAcc;
658 fCountRecoITSClusters+= icountRecoITSClusters;
659 fCountRecoPPR+= icountRecoPPR;
660 fCountRecoPID+= icountRecoPID;
661
662 fHistEventsProcessed->Fill(0);
f2703bd2 663
664 delete[] containerInput;
665 containerInput = 0x0;
666 delete[] containerInputMC;
667 containerInputMC = 0x0;
668 delete cfVtxHF;
669
379592af 670}
671
672//___________________________________________________________________________
673void AliCFTaskVertexingHF::Terminate(Option_t*)
674{
675 // The Terminate() function is the last function to be called during
676 // a query. It always runs on the client, it can be used to present
677 // the results graphically or save the results to file.
678
679 AliAnalysisTaskSE::Terminate();
680
681 AliInfo(Form("Found %i MC particles that are D0 in MC, in %d events",fCountMC,fEvents));
682 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, in %d events",fCountAcc,fEvents));
683 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fEvents));
684 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, and satisfy ITS+TPC refit requirementin %d events",fCountRefit,fEvents));
685 AliInfo(Form("Found %i reco D0 that are decaying in K+pi, in %d events",fCountReco,fEvents));
686 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and are in the requested acceptance, in %d events",fCountRecoAcc,fEvents));
687 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and have at least 5 clusters in ITS, in %d events",fCountRecoITSClusters,fEvents));
688 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and satisfy PPR cuts, in %d events",fCountRecoPPR,fEvents));
689
690 // draw some example plots....
691
692 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
693 if(!cont) {
694 printf("CONTAINER NOT FOUND\n");
695 return;
696 }
697 // projecting the containers to obtain histograms
698 // first argument = variable, second argument = step
699
700 // MC-level
701 TH1D* h00 = cont->ShowProjection(0,0) ; // pt
702 TH1D* h10 = cont->ShowProjection(1,0) ; // rapidity
703 TH1D* h20 = cont->ShowProjection(2,0) ; // cosThetaStar
704 TH1D* h30 = cont->ShowProjection(3,0) ; // pTpi
705 TH1D* h40 = cont->ShowProjection(4,0) ; // pTK
706 TH1D* h50 = cont->ShowProjection(5,0) ; // cT
707 TH1D* h60 = cont->ShowProjection(6,0) ; // dca
708 TH1D* h70 = cont->ShowProjection(7,0) ; // d0pi
709 TH1D* h80 = cont->ShowProjection(8,0) ; // d0K
710 TH1D* h90 = cont->ShowProjection(9,0) ; // d0xd0
711 TH1D* h100 = cont->ShowProjection(10,0) ; // cosPointingAngle
712 TH1D* h110 = cont->ShowProjection(11,0) ; // phi
713
714 // MC-Acceptance level
715 TH1D* h01 = cont->ShowProjection(0,1) ; // pt
716 TH1D* h11 = cont->ShowProjection(1,1) ; // rapidity
717 TH1D* h21 = cont->ShowProjection(2,1) ; // cosThetaStar
718 TH1D* h31 = cont->ShowProjection(3,1) ; // pTpi
719 TH1D* h41 = cont->ShowProjection(4,1) ; // pTK
720 TH1D* h51 = cont->ShowProjection(5,1) ; // cT
721 TH1D* h61 = cont->ShowProjection(6,1) ; // dca
722 TH1D* h71 = cont->ShowProjection(7,1) ; // d0pi
723 TH1D* h81 = cont->ShowProjection(8,1) ; // d0K
724 TH1D* h91 = cont->ShowProjection(9,1) ; // d0xd0
725 TH1D* h101 = cont->ShowProjection(10,1) ; // cosPointingAngle
726 TH1D* h111 = cont->ShowProjection(11,1) ; // phi
727
728 // Reco-level
729 TH1D* h04 = cont->ShowProjection(0,4) ; // pt
730 TH1D* h14 = cont->ShowProjection(1,4) ; // rapidity
731 TH1D* h24 = cont->ShowProjection(2,4) ; // cosThetaStar
732 TH1D* h34 = cont->ShowProjection(3,4) ; // pTpi
733 TH1D* h44 = cont->ShowProjection(4,4) ; // pTK
734 TH1D* h54 = cont->ShowProjection(5,4) ; // cT
735 TH1D* h64 = cont->ShowProjection(6,4) ; // dca
736 TH1D* h74 = cont->ShowProjection(7,4) ; // d0pi
737 TH1D* h84 = cont->ShowProjection(8,4) ; // d0K
738 TH1D* h94 = cont->ShowProjection(9,4) ; // d0xd0
739 TH1D* h104 = cont->ShowProjection(10,4) ; // cosPointingAngle
740 TH1D* h114 = cont->ShowProjection(11,4) ; // phi
741
742 h00->SetTitle("pT_D0 (GeV/c)");
743 h10->SetTitle("rapidity");
744 h20->SetTitle("cosThetaStar");
745 h30->SetTitle("pT_pi (GeV/c)");
746 h40->SetTitle("pT_K (Gev/c)");
747 h50->SetTitle("cT (#mum)");
748 h60->SetTitle("dca (#mum)");
749 h70->SetTitle("d0_pi (#mum)");
750 h80->SetTitle("d0_K (#mum)");
751 h90->SetTitle("d0xd0 (#mum^2)");
752 h100->SetTitle("cosPointingAngle");
753 h100->SetTitle("phi (rad)");
754
755 h00->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
756 h10->GetXaxis()->SetTitle("rapidity");
757 h20->GetXaxis()->SetTitle("cosThetaStar");
758 h30->GetXaxis()->SetTitle("pT_pi (GeV/c)");
759 h40->GetXaxis()->SetTitle("pT_K (Gev/c)");
760 h50->GetXaxis()->SetTitle("cT (#mum)");
761 h60->GetXaxis()->SetTitle("dca (#mum)");
762 h70->GetXaxis()->SetTitle("d0_pi (#mum)");
763 h80->GetXaxis()->SetTitle("d0_K (#mum)");
764 h90->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
765 h100->GetXaxis()->SetTitle("cosPointingAngle");
766 h110->GetXaxis()->SetTitle("phi (rad)");
767
768 h01->SetTitle("pT_D0 (GeV/c)");
769 h11->SetTitle("rapidity");
770 h21->SetTitle("cosThetaStar");
771 h31->SetTitle("pT_pi (GeV/c)");
772 h41->SetTitle("pT_K (Gev/c)");
773 h51->SetTitle("cT (#mum)");
774 h61->SetTitle("dca (#mum)");
775 h71->SetTitle("d0_pi (#mum)");
776 h81->SetTitle("d0_K (#mum)");
777 h91->SetTitle("d0xd0 (#mum^2)");
778 h101->SetTitle("cosPointingAngle");
779 h111->GetXaxis()->SetTitle("phi (rad)");
780
781 h01->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
782 h11->GetXaxis()->SetTitle("rapidity");
783 h21->GetXaxis()->SetTitle("cosThetaStar");
784 h31->GetXaxis()->SetTitle("pT_pi (GeV/c)");
785 h41->GetXaxis()->SetTitle("pT_K (Gev/c)");
786 h51->GetXaxis()->SetTitle("cT (#mum)");
787 h61->GetXaxis()->SetTitle("dca (#mum)");
788 h71->GetXaxis()->SetTitle("d0_pi (#mum)");
789 h81->GetXaxis()->SetTitle("d0_K (#mum)");
790 h91->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
791 h101->GetXaxis()->SetTitle("cosPointingAngle");
792 h111->GetXaxis()->SetTitle("phi (rad)");
793
794 h04->SetTitle("pT_D0 (GeV/c)");
795 h14->SetTitle("rapidity");
796 h24->SetTitle("cosThetaStar");
797 h34->SetTitle("pT_pi (GeV/c)");
798 h44->SetTitle("pT_K (Gev/c)");
799 h54->SetTitle("cT (#mum)");
800 h64->SetTitle("dca (#mum)");
801 h74->SetTitle("d0_pi (#mum)");
802 h84->SetTitle("d0_K (#mum)");
803 h94->SetTitle("d0xd0 (#mum^2)");
804 h104->SetTitle("cosPointingAngle");
805 h114->GetXaxis()->SetTitle("phi (rad)");
806
807 h04->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
808 h14->GetXaxis()->SetTitle("rapidity");
809 h24->GetXaxis()->SetTitle("cosThetaStar");
810 h34->GetXaxis()->SetTitle("pT_pi (GeV/c)");
811 h44->GetXaxis()->SetTitle("pT_K (Gev/c)");
812 h54->GetXaxis()->SetTitle("cT (#mum)");
813 h64->GetXaxis()->SetTitle("dca (#mum)");
814 h74->GetXaxis()->SetTitle("d0_pi (#mum)");
815 h84->GetXaxis()->SetTitle("d0_K (#mum)");
816 h94->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
817 h104->GetXaxis()->SetTitle("cosPointingAngle");
818 h114->GetXaxis()->SetTitle("phi (rad)");
819
820 Double_t max0 = h00->GetMaximum();
821 Double_t max1 = h10->GetMaximum();
822 Double_t max2 = h20->GetMaximum();
823 Double_t max3 = h30->GetMaximum();
824 Double_t max4 = h40->GetMaximum();
825 Double_t max5 = h50->GetMaximum();
826 Double_t max6 = h60->GetMaximum();
827 Double_t max7 = h70->GetMaximum();
828 Double_t max8 = h80->GetMaximum();
829 Double_t max9 = h90->GetMaximum();
830 Double_t max10 = h100->GetMaximum();
831 Double_t max11 = h110->GetMaximum();
832
833 h00->GetYaxis()->SetRangeUser(0,max0*1.2);
834 h10->GetYaxis()->SetRangeUser(0,max1*1.2);
835 h20->GetYaxis()->SetRangeUser(0,max2*1.2);
836 h30->GetYaxis()->SetRangeUser(0,max3*1.2);
837 h40->GetYaxis()->SetRangeUser(0,max4*1.2);
838 h50->GetYaxis()->SetRangeUser(0,max5*1.2);
839 h60->GetYaxis()->SetRangeUser(0,max6*1.2);
840 h70->GetYaxis()->SetRangeUser(0,max7*1.2);
841 h80->GetYaxis()->SetRangeUser(0,max8*1.2);
842 h90->GetYaxis()->SetRangeUser(0,max9*1.2);
843 h100->GetYaxis()->SetRangeUser(0,max10*1.2);
844 h110->GetYaxis()->SetRangeUser(0,max11*1.2);
845
846 h01->GetYaxis()->SetRangeUser(0,max0*1.2);
847 h11->GetYaxis()->SetRangeUser(0,max1*1.2);
848 h21->GetYaxis()->SetRangeUser(0,max2*1.2);
849 h31->GetYaxis()->SetRangeUser(0,max3*1.2);
850 h41->GetYaxis()->SetRangeUser(0,max4*1.2);
851 h51->GetYaxis()->SetRangeUser(0,max5*1.2);
852 h61->GetYaxis()->SetRangeUser(0,max6*1.2);
853 h71->GetYaxis()->SetRangeUser(0,max7*1.2);
854 h81->GetYaxis()->SetRangeUser(0,max8*1.2);
855 h91->GetYaxis()->SetRangeUser(0,max9*1.2);
856 h101->GetYaxis()->SetRangeUser(0,max10*1.2);
857 h111->GetYaxis()->SetRangeUser(0,max11*1.2);
858
859 h00->SetMarkerStyle(20);
860 h10->SetMarkerStyle(24);
861 h20->SetMarkerStyle(21);
862 h30->SetMarkerStyle(25);
863 h40->SetMarkerStyle(27);
864 h50->SetMarkerStyle(28);
865 h60->SetMarkerStyle(20);
866 h70->SetMarkerStyle(24);
867 h80->SetMarkerStyle(21);
868 h90->SetMarkerStyle(25);
869 h100->SetMarkerStyle(27);
870 h110->SetMarkerStyle(28);
871
872 h00->SetMarkerColor(2);
873 h10->SetMarkerColor(2);
874 h20->SetMarkerColor(2);
875 h30->SetMarkerColor(2);
876 h40->SetMarkerColor(2);
877 h50->SetMarkerColor(2);
878 h60->SetMarkerColor(2);
879 h70->SetMarkerColor(2);
880 h80->SetMarkerColor(2);
881 h90->SetMarkerColor(2);
882 h100->SetMarkerColor(2);
883 h110->SetMarkerColor(2);
884
885 h01->SetMarkerStyle(20) ;
886 h11->SetMarkerStyle(24) ;
887 h21->SetMarkerStyle(21) ;
888 h31->SetMarkerStyle(25) ;
889 h41->SetMarkerStyle(27) ;
890 h51->SetMarkerStyle(28) ;
891 h61->SetMarkerStyle(20);
892 h71->SetMarkerStyle(24);
893 h81->SetMarkerStyle(21);
894 h91->SetMarkerStyle(25);
895 h101->SetMarkerStyle(27);
896 h111->SetMarkerStyle(28);
897
898 h01->SetMarkerColor(8);
899 h11->SetMarkerColor(8);
900 h21->SetMarkerColor(8);
901 h31->SetMarkerColor(8);
902 h41->SetMarkerColor(8);
903 h51->SetMarkerColor(8);
904 h61->SetMarkerColor(8);
905 h71->SetMarkerColor(8);
906 h81->SetMarkerColor(8);
907 h91->SetMarkerColor(8);
908 h101->SetMarkerColor(8);
909 h111->SetMarkerColor(8);
910
911 h04->SetMarkerStyle(20) ;
912 h14->SetMarkerStyle(24) ;
913 h24->SetMarkerStyle(21) ;
914 h34->SetMarkerStyle(25) ;
915 h44->SetMarkerStyle(27) ;
916 h54->SetMarkerStyle(28) ;
917 h64->SetMarkerStyle(20);
918 h74->SetMarkerStyle(24);
919 h84->SetMarkerStyle(21);
920 h94->SetMarkerStyle(25);
921 h104->SetMarkerStyle(27);
922 h114->SetMarkerStyle(28);
923
924 h04->SetMarkerColor(4);
925 h14->SetMarkerColor(4);
926 h24->SetMarkerColor(4);
927 h34->SetMarkerColor(4);
928 h44->SetMarkerColor(4);
929 h54->SetMarkerColor(4);
930 h64->SetMarkerColor(4);
931 h74->SetMarkerColor(4);
932 h84->SetMarkerColor(4);
933 h94->SetMarkerColor(4);
934 h104->SetMarkerColor(4);
935 h114->SetMarkerColor(4);
936
937 gStyle->SetCanvasColor(0);
938 gStyle->SetFrameFillColor(0);
939 gStyle->SetTitleFillColor(0);
940 gStyle->SetStatColor(0);
941
942 // drawing in 2 separate canvas for a matter of clearity
943 TCanvas * c1 =new TCanvas("c1New","pT, rapidiy, cosThetaStar",1100,1600);
944 c1->Divide(3,3);
945
946 c1->cd(1);
947 h00->Draw("p");
948 c1->cd(1);
949 c1->cd(2);
950 h01->Draw("p");
951 c1->cd(2);
952 c1->cd(3);
953 h04->Draw("p");
954 c1->cd(3);
955 c1->cd(4);
956 h10->Draw("p");
957 c1->cd(4);
958 c1->cd(5);
959 h11->Draw("p");
960 c1->cd(5);
961 c1->cd(6);
962 h14->Draw("p");
963 c1->cd(6);
964 c1->cd(7);
965 h20->Draw("p");
966 c1->cd(7);
967 c1->cd(8);
968 h21->Draw("p");
969 c1->cd(8);
970 c1->cd(9);
971 h24->Draw("p");
972 c1->cd(9);
973 c1->cd();
974
975 TCanvas * c2 =new TCanvas("c2New","pTpi, pTK, cT",1100,1600);
976 c2->Divide(3,3);
977 c2->cd(1);
978 h30->Draw("p");
979 c2->cd(1);
980 c2->cd(2);
981 h31->Draw("p");
982 c2->cd(2);
983 c2->cd(3);
984 h34->Draw("p");
985 c2->cd(3);
986 c2->cd(4);
987 h40->Draw("p");
988 c2->cd(4);
989 c2->cd(5);
990 h41->Draw("p");
991 c2->cd(5);
992 c2->cd(6);
993 h44->Draw("p");
994 c2->cd(6);
995 c2->cd(7);
996 h50->Draw("p");
997 c2->cd(7);
998 c2->cd(8);
999 h51->Draw("p");
1000 c2->cd(8);
1001 c2->cd(9);
1002 h54->Draw("p");
1003 c2->cd(9);
1004 c2->cd();
1005
1006 TCanvas * c3 =new TCanvas("c3New","dca, d0pi, d0K",1100,1600);
1007 c3->Divide(3,3);
1008 c3->cd(1);
1009 h60->Draw("p");
1010 c3->cd(1);
1011 c3->cd(2);
1012 h61->Draw("p");
1013 c3->cd(2);
1014 c3->cd(3);
1015 h64->Draw("p");
1016 c3->cd(3);
1017 c3->cd(4);
1018 h70->Draw("p");
1019 c3->cd(4);
1020 c3->cd(5);
1021 h71->Draw("p");
1022 c3->cd(5);
1023 c3->cd(6);
1024 h74->Draw("p");
1025 c3->cd(6);
1026 c3->cd(7);
1027 h80->Draw("p");
1028 c3->cd(7);
1029 c3->cd(8);
1030 h81->Draw("p");
1031 c3->cd(8);
1032 c3->cd(9);
1033 h84->Draw("p");
1034 c3->cd(9);
1035 c3->cd();
1036
1037 TCanvas * c4 =new TCanvas("c4New","d0xd0, cosPointingAngle, phi",1100,1600);
1038 c4->Divide(3,3);
1039 c4->cd(1);
1040 h90->Draw("p");
1041 c4->cd(1);
1042 c4->cd(2);
1043 h91->Draw("p");
1044 c4->cd(2);
1045 c4->cd(3);
1046 h94->Draw("p");
1047 c4->cd(3);
1048 c4->cd(4);
1049 h100->Draw("p");
1050 c4->cd(4);
1051 c4->cd(5);
1052 h101->Draw("p");
1053 c4->cd(5);
1054 c4->cd(6);
1055 h104->Draw("p");
1056 c4->cd(6);
1057 c4->cd(7);
1058 h110->Draw("p");
1059 c4->cd(7);
1060 c4->cd(8);
1061 h111->Draw("p");
1062 c4->cd(8);
1063 c4->cd(9);
1064 h114->Draw("p");
1065 c4->cd(9);
1066 c4->cd();
1067
1068 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
1069
1070 TH2D* corr1 =hcorr->Projection(0,2);
1071 TH2D* corr2 = hcorr->Projection(1,3);
1072
1073 TCanvas * c7 =new TCanvas("c7New","",800,400);
1074 c7->Divide(2,1);
1075 c7->cd(1);
1076 corr1->Draw("text");
1077 c7->cd(2);
1078 corr2->Draw("text");
1079
1080
1081 TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
1082
1083 corr1->Write();
1084 corr2->Write();
1085 h00->Write("pT_D0_step0");
1086 h10->Write("rapidity_step0");
1087 h20->Write("cosThetaStar_step0");
1088 h30->Write("pT_pi_step0");
1089 h40->Write("pT_K_step0");
1090 h50->Write("cT_step0");
1091 h60->Write("dca_step0");
1092 h70->Write("d0_pi_step0");
1093 h80->Write("d0_K_step0");
1094 h90->Write("d0xd0_step0");
1095 h100->Write("cosPointingAngle_step0");
1096 h110->Write("phi_step0");
1097
1098 h01->Write("pT_D0_step1");
1099 h11->Write("rapidity_step1");
1100 h21->Write("cosThetaStar_step1");
1101 h31->Write("pT_pi_step1");
1102 h41->Write("pT_K_step1");
1103 h51->Write("cT_step1");
1104 h61->Write("dca_step1");
1105 h71->Write("d0_pi_step1");
1106 h81->Write("d0_K_step1");
1107 h91->Write("d0xd0_step1");
1108 h101->Write("cosPointingAngle_step1");
1109 h111->Write("phi_step1");
1110
1111 h04->Write("pT_D0_step2");
1112 h14->Write("rapidity_step2");
1113 h24->Write("cosThetaStar_step2");
1114 h34->Write("pT_pi_step2");
1115 h44->Write("pT_K_step2");
1116 h54->Write("cT_step2");
1117 h64->Write("dca_step2");
1118 h74->Write("d0_pi_step2");
1119 h80->Write("d0_K_step2");
1120 h94->Write("d0xd0_step2");
1121 h104->Write("cosPointingAngle_step2");
1122 h114->Write("phi_step2");
1123
1124 file_projection->Close();
1125
1126
1127
1128 /*
1129 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
1130 c2->SaveAs("Plots/pTpi_pTK_cT.eps");
1131 c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
1132 c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
1133
1134 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
1135 c2->SaveAs("Plots/pTpi_pTK_cT.gif");
1136 c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
1137 c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
1138 */
1139}
1140
1141//___________________________________________________________________________
f2703bd2 1142void AliCFTaskVertexingHF::UserCreateOutputObjects()
1143{
379592af 1144 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1145 //TO BE SET BEFORE THE EXECUTION OF THE TASK
1146 //
1147 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1148
1149 //slot #1
1150 OpenFile(1);
1151 fHistEventsProcessed = new TH1I("CFHFchist0","",1,0,1) ;
1152}
1153
f2703bd2 1154
1155//_________________________________________________________________________
1156Double_t AliCFTaskVertexingHF::GetWeight(Float_t pt)
1157{
1158 //
1159 // calculating the weight to fill the container
1160 //
1161
1162 // FNOLL central:
1163 // p0 = 1.63297e-01 --> 0.322643
1164 // p1 = 2.96275e+00
1165 // p2 = 2.30301e+00
1166 // p3 = 2.50000e+00
1167
1168 // PYTHIA
1169 // p0 = 1.85906e-01 --> 0.36609
1170 // p1 = 1.94635e+00
1171 // p2 = 1.40463e+00
1172 // p3 = 2.50000e+00
1173
1174 Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
1175 Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
1176
1177 Double_t dndpt_func1 = dNdptFit(pt,func1);
1178 Double_t dndpt_func2 = dNdptFit(pt,func2);
1179 AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
1180 return dndpt_func1/dndpt_func2;
1181}
1182
1183//__________________________________________________________________________________________________
1184Double_t AliCFTaskVertexingHF::dNdptFit(Float_t pt, Double_t* par)
1185{
1186 //
1187 // calculating dNdpt
1188 //
1189
1190 Double_t denom = TMath::Power((pt/par[1]), par[3] );
1191 Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
1192
1193 return dNdpt;
1194}