1 /**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
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 //-----------------------------------------------------------------------
33 #include <TParticle.h>
34 #include <TDatabasePDG.h>
39 #include "AliCFTaskVertexingHF.h"
41 #include "AliMCEvent.h"
42 #include "AliCFManager.h"
43 #include "AliCFContainer.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"
51 #include "AliAODRecoDecayHF3Prong.h"
52 #include "AliAODRecoDecayHF4Prong.h"
53 #include "AliAODRecoCascadeHF.h"
54 #include "AliAODMCParticle.h"
55 #include "AliAODMCHeader.h"
56 #include "AliESDtrack.h"
58 #include "THnSparse.h"
60 #include "AliESDtrackCuts.h"
61 #include "AliRDHFCuts.h"
62 #include "AliRDHFCutsD0toKpi.h"
63 #include "AliRDHFCutsDplustoKpipi.h"
64 #include "AliRDHFCutsDStartoKpipi.h"
65 #include "AliRDHFCutsDstoKKpi.h"
66 #include "AliRDHFCutsLctopKpi.h"
67 #include "AliRDHFCutsD0toKpipipi.h"
68 #include "AliCFVertexingHF2Prong.h"
69 #include "AliCFVertexingHF3Prong.h"
70 #include "AliCFVertexingHF.h"
71 #include "AliAnalysisDataSlot.h"
72 #include "AliAnalysisDataContainer.h"
74 //__________________________________________________________________________
75 AliCFTaskVertexingHF::AliCFTaskVertexingHF() :
78 fHistEventsProcessed(0x0),
86 fCountRecoITSClusters(0),
91 fFillFromGenerated(kFALSE),
93 fAcceptanceUnf(kTRUE),
106 //___________________________________________________________________________
107 AliCFTaskVertexingHF::AliCFTaskVertexingHF(const Char_t* name, AliRDHFCuts* cuts) :
108 AliAnalysisTaskSE(name),
110 fHistEventsProcessed(0x0),
118 fCountRecoITSClusters(0),
123 fFillFromGenerated(kFALSE),
124 fOriginDselection(0),
125 fAcceptanceUnf(kTRUE),
135 // Constructor. Initialization of Inputs and Outputs
138 DefineInput(0) and DefineOutput(0)
139 are taken care of by AliAnalysisTaskSE constructor
141 DefineOutput(1,TH1I::Class());
142 DefineOutput(2,AliCFContainer::Class());
143 DefineOutput(3,THnSparseD::Class());
144 DefineOutput(4,AliRDHFCuts::Class());
149 //___________________________________________________________________________
150 AliCFTaskVertexingHF& AliCFTaskVertexingHF::operator=(const AliCFTaskVertexingHF& c)
153 // Assignment operator
156 AliAnalysisTaskSE::operator=(c) ;
157 fCFManager = c.fCFManager;
158 fHistEventsProcessed = c.fHistEventsProcessed;
164 //___________________________________________________________________________
165 AliCFTaskVertexingHF::AliCFTaskVertexingHF(const AliCFTaskVertexingHF& c) :
166 AliAnalysisTaskSE(c),
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),
180 fDecayChannel(c.fDecayChannel),
181 fFillFromGenerated(c.fFillFromGenerated),
182 fOriginDselection(c.fOriginDselection),
183 fAcceptanceUnf(c.fAcceptanceUnf),
185 fUseWeight(c.fUseWeight),
188 fPartName(c.fPartName),
189 fDauNames(c.fDauNames),
197 //___________________________________________________________________________
198 AliCFTaskVertexingHF::~AliCFTaskVertexingHF()
203 if (fCFManager) delete fCFManager ;
204 if (fHistEventsProcessed) delete fHistEventsProcessed ;
205 if (fCorrelation) delete fCorrelation ;
206 if (fCuts) delete fCuts;
209 //_________________________________________________________________________-
210 void AliCFTaskVertexingHF::Init()
216 if (fDebug>1) printf("AliCFTaskVertexingHF::Init()");
217 AliRDHFCuts *copyfCuts = 0x0;
219 switch (fDecayChannel){
221 copyfCuts = new AliRDHFCutsD0toKpi(*(dynamic_cast<AliRDHFCutsD0toKpi*>(fCuts)));
228 copyfCuts = new AliRDHFCutsDStartoKpipi(*(dynamic_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
235 copyfCuts = new AliRDHFCutsDplustoKpipi(*(dynamic_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
242 copyfCuts = new AliRDHFCutsLctopKpi(*(dynamic_cast<AliRDHFCutsLctopKpi*>(fCuts)));
249 copyfCuts = new AliRDHFCutsDstoKKpi(*(dynamic_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
256 copyfCuts = new AliRDHFCutsD0toKpipipi(*(dynamic_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
259 fDauNames="K+pi+pi+pi";
263 AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exiting...");
267 const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
268 copyfCuts->SetName(nameoutput);
271 PostData(4, copyfCuts);
276 //_________________________________________________
277 void AliCFTaskVertexingHF::UserExec(Option_t *)
280 // Main loop function
283 PostData(1,fHistEventsProcessed) ;
284 PostData(2,fCFManager->GetParticleContainer()) ;
285 PostData(3,fCorrelation) ;
287 AliESDtrackCuts* trackCuts = fCuts->GetTrackCuts();
289 if (fFillFromGenerated){
290 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
294 Error("UserExec","NO EVENT FOUND!");
298 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
300 TClonesArray *arrayBranch=0;
302 if(!aodEvent && AODEvent() && IsStandardAOD()) {
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();
314 switch (fDecayChannel){
316 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
320 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
326 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
330 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
339 switch (fDecayChannel){
341 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
345 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
351 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
355 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm4Prong");
363 AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
367 AliError("Could not find array of HF vertices");
373 fCFManager->SetRecEventInfo(aodEvent);
374 fCFManager->SetMCEventInfo(aodEvent);
376 //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro.
378 Double_t* containerInput = new Double_t[fNvar];
379 Double_t* containerInputMC = new Double_t[fNvar];
381 //loop on the MC event
383 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
385 AliError("Could not find Monte-Carlo in AOD");
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;
399 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
401 AliError("Could not find MC Header in AOD");
405 AliCFVertexingHF* cfVtxHF=0x0;
406 switch (fDecayChannel){
408 cfVtxHF = new AliCFVertexingHF2Prong(mcArray, fOriginDselection);
412 // cfVtxHF = new AliCFVertexingHFCascade(mcArray, originDselection); // not there yet
418 cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel);
422 //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection); // not there yet
429 Double_t zPrimVertex = aodVtx ->GetZ();
430 Double_t zMCVertex = mcHeader->GetVtxZ();
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);
437 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
439 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
441 // check the MC-level cuts, must be the desidered particle
442 if (!fCFManager->CheckParticleCuts(0, mcPart)) continue; // 0 stands for MC level
444 cfVtxHF->SetMCCandidateParam(iPart);
445 cfVtxHF->SetNVar(fNvar);
447 cquarks += cfVtxHF->MCcquarkCounting(mcPart);
449 if (!(cfVtxHF->SetLabelArray())){
450 AliDebug(2,Form("Impossible to set the label array (decaychannel = %d)",fDecayChannel));
454 //check the candiate family at MC level
455 if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) {
456 AliDebug(2,Form("Check on the family wrong!!! (decaychannel = %d)",fDecayChannel));
460 AliDebug(2,Form("Check on the family OK!!! (decaychannel = %d)",fDecayChannel));
463 //Fill the MC container
464 Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
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");
475 fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight);
477 AliDebug(3,"MC cointainer filled \n");
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();
484 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight);
485 AliDebug(3,"MC acceptance cut passed\n");
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");
495 //mc Refit requirement
496 Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, trackCuts);
498 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight);
499 AliDebug(3,"MC Refit cut passed and container filled\n");
503 AliDebug(3,"MC Refit cut not passed\n");
508 AliDebug (3, "MC vertex step not passed\n");
513 AliDebug (3, "MC in acceptance step not passed\n");
518 AliDebug (3, "MC container not filled\n");
522 if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks));
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()));
528 // Now go to rec level
529 fCountMC += icountMC;
530 fCountAcc += icountAcc;
531 fCountVertex+= icountVertex;
532 fCountRefit+= icountRefit;
534 AliDebug(2,Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel));
536 for(Int_t iCandid = 0; iCandid<arrayBranch->GetEntriesFast();iCandid++){
537 AliAODRecoDecayHF* charmCandidate=0x0;
538 switch (fDecayChannel){
540 charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid);
544 charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid);
550 charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid);
554 charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid);
561 Bool_t unsetvtx=kFALSE;
562 if(!charmCandidate->GetOwnPrimaryVtx()) {
563 charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
567 Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate);
568 if (!signAssociation){
569 charmCandidate = 0x0;
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));
579 Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
582 if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue;
585 Bool_t recoStep = cfVtxHF->RecoStep();
586 Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
588 if (recoStep && recoContFilled && vtxCheck){
589 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
591 AliDebug(3,"Reco step passed and container filled\n");
593 //Reco in the acceptance -- take care of UNFOLDING!!!!
594 Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(trackCuts);
595 if (recoAcceptanceStep) {
596 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
598 AliDebug(3,"Reco acceptance cut passed and container filled\n");
601 Double_t fill[4]; //fill response matrix
602 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
603 if (bUnfolding) fCorrelation->Fill(fill);
606 //Number of ITS cluster requirements
607 Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
608 if (recoITSnCluster){
609 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
610 icountRecoITSClusters++;
611 AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
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);
618 if (recoAnalysisCuts > 3){ // Ds case, where more possibilities are considered
619 if (recoAnalysisCuts >= 8){
620 recoAnalysisCuts -= 8; // removing K0star mass
622 if (recoAnalysisCuts >= 4){
623 recoAnalysisCuts -= 4; // removing Phi mass
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);
631 AliDebug(3,"Reco Analysis cuts passed and container filled \n");
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);
639 AliDebug(3,"Reco PID cuts passed and container filled \n");
641 Double_t fill[4]; //fill response matrix
642 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
643 if (bUnfolding) fCorrelation->Fill(fill);
647 AliDebug(3, "Analysis Cuts step not passed \n");
652 AliDebug(3, "PID selection not passed \n");
657 AliDebug(3, "Number of ITS cluster step not passed\n");
662 AliDebug(3, "Reco acceptance step not passed\n");
667 AliDebug(3, "Reco step not passed\n");
672 if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
673 } // end loop on candidate
675 fCountReco+= icountReco;
676 fCountRecoAcc+= icountRecoAcc;
677 fCountRecoITSClusters+= icountRecoITSClusters;
678 fCountRecoPPR+= icountRecoPPR;
679 fCountRecoPID+= icountRecoPID;
681 fHistEventsProcessed->Fill(0);
683 delete[] containerInput;
684 containerInput = 0x0;
685 delete[] containerInputMC;
686 containerInputMC = 0x0;
691 //___________________________________________________________________________
692 void AliCFTaskVertexingHF::Terminate(Option_t*)
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.
698 AliAnalysisTaskSE::Terminate();
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));
710 // draw some example plots....
712 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
714 printf("CONTAINER NOT FOUND\n");
717 // projecting the containers to obtain histograms
718 // first argument = variable, second argument = step
721 for(Int_t iC=0;iC<12; iC++){
723 h[0][iC] = cont->ShowProjection(iC,0);
724 // MC-Acceptance level
725 h[1][iC] = cont->ShowProjection(iC,1);
727 h[2][iC] = cont->ShowProjection(iC,4);
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)";
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)";
757 Int_t markers[12]={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]);
773 gStyle->SetCanvasColor(0);
774 gStyle->SetFrameFillColor(0);
775 gStyle->SetTitleFillColor(0);
776 gStyle->SetStatColor(0);
778 // drawing in 2 separate canvas for a matter of clearity
779 TCanvas * c1 =new TCanvas("c1New","Vars 0,1,2",1100,1600);
782 for(Int_t iVar=0; iVar<3; iVar++){
784 h[0][iVar]->Draw("p");
786 h[1][iVar]->Draw("p");
788 h[2][iVar]->Draw("p");
791 TCanvas * c2 =new TCanvas("c2New","Vars 3,4,5",1100,1600);
794 for(Int_t iVar=3; iVar<6; iVar++){
796 h[0][iVar]->Draw("p");
798 h[1][iVar]->Draw("p");
800 h[2][iVar]->Draw("p");
804 TCanvas * c3 =new TCanvas("c3New","Vars 6,7,8",1100,1600);
807 for(Int_t iVar=6; iVar<9; iVar++){
809 h[0][iVar]->Draw("p");
811 h[1][iVar]->Draw("p");
813 h[2][iVar]->Draw("p");
817 TCanvas * c4 =new TCanvas("c4New","Vars 9,10,11",1100,1600);
820 for(Int_t iVar=9; iVar<11; iVar++){
822 h[0][iVar]->Draw("p");
824 h[1][iVar]->Draw("p");
826 h[2][iVar]->Draw("p");
830 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
832 TH2D* corr1 =hcorr->Projection(0,2);
833 TH2D* corr2 = hcorr->Projection(1,3);
835 TCanvas * c7 =new TCanvas("c7New","",800,400);
843 TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
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()));
852 file_projection->Close();
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");
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");
869 //___________________________________________________________________________
870 void AliCFTaskVertexingHF::UserCreateOutputObjects()
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
875 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
879 fHistEventsProcessed = new TH1I("CFHFchist0","",1,0,1) ;
883 //_________________________________________________________________________
884 Double_t AliCFTaskVertexingHF::GetWeight(Float_t pt)
887 // calculating the weight to fill the container
891 // p0 = 1.63297e-01 --> 0.322643
897 // p0 = 1.85906e-01 --> 0.36609
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};
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;
911 //__________________________________________________________________________________________________
912 Double_t AliCFTaskVertexingHF::dNdptFit(Float_t pt, Double_t* par)
918 Double_t denom = TMath::Power((pt/par[1]), par[3] );
919 Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);