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 "AliCFVertexingHFCascade.h"
71 #include "AliCFVertexingHF.h"
72 #include "AliAnalysisDataSlot.h"
73 #include "AliAnalysisDataContainer.h"
75 //__________________________________________________________________________
76 AliCFTaskVertexingHF::AliCFTaskVertexingHF() :
79 fHistEventsProcessed(0x0),
87 fCountRecoITSClusters(0),
92 fFillFromGenerated(kFALSE),
94 fAcceptanceUnf(kTRUE),
102 fCentralitySelection(kTRUE),
104 fRejectIfNoQuark(kTRUE),
105 fUseMCVertex(kFALSE),
112 //___________________________________________________________________________
113 AliCFTaskVertexingHF::AliCFTaskVertexingHF(const Char_t* name, AliRDHFCuts* cuts) :
114 AliAnalysisTaskSE(name),
116 fHistEventsProcessed(0x0),
124 fCountRecoITSClusters(0),
129 fFillFromGenerated(kFALSE),
130 fOriginDselection(0),
131 fAcceptanceUnf(kTRUE),
139 fCentralitySelection(kTRUE),
141 fRejectIfNoQuark(kTRUE),
142 fUseMCVertex(kFALSE),
146 // Constructor. Initialization of Inputs and Outputs
149 DefineInput(0) and DefineOutput(0)
150 are taken care of by AliAnalysisTaskSE constructor
152 DefineOutput(1,TH1I::Class());
153 DefineOutput(2,AliCFContainer::Class());
154 DefineOutput(3,THnSparseD::Class());
155 DefineOutput(4,AliRDHFCuts::Class());
160 //___________________________________________________________________________
161 AliCFTaskVertexingHF& AliCFTaskVertexingHF::operator=(const AliCFTaskVertexingHF& c)
164 // Assignment operator
167 AliAnalysisTaskSE::operator=(c) ;
168 fCFManager = c.fCFManager;
169 fHistEventsProcessed = c.fHistEventsProcessed;
175 //___________________________________________________________________________
176 AliCFTaskVertexingHF::AliCFTaskVertexingHF(const AliCFTaskVertexingHF& c) :
177 AliAnalysisTaskSE(c),
178 fCFManager(c.fCFManager),
179 fHistEventsProcessed(c.fHistEventsProcessed),
180 fCorrelation(c.fCorrelation),
181 fCountMC(c.fCountMC),
182 fCountAcc(c.fCountAcc),
183 fCountVertex(c.fCountVertex),
184 fCountRefit(c.fCountRefit),
185 fCountReco(c.fCountReco),
186 fCountRecoAcc(c.fCountRecoAcc),
187 fCountRecoITSClusters(c.fCountRecoITSClusters),
188 fCountRecoPPR(c.fCountRecoPPR),
189 fCountRecoPID(c.fCountRecoPID),
191 fDecayChannel(c.fDecayChannel),
192 fFillFromGenerated(c.fFillFromGenerated),
193 fOriginDselection(c.fOriginDselection),
194 fAcceptanceUnf(c.fAcceptanceUnf),
196 fUseWeight(c.fUseWeight),
199 fPartName(c.fPartName),
200 fDauNames(c.fDauNames),
202 fCentralitySelection(c.fCentralitySelection),
203 fFakeSelection(c.fFakeSelection),
204 fRejectIfNoQuark(c.fRejectIfNoQuark),
205 fUseMCVertex(c.fUseMCVertex),
206 fDsOption(c.fDsOption)
213 //___________________________________________________________________________
214 AliCFTaskVertexingHF::~AliCFTaskVertexingHF()
219 if (fCFManager) delete fCFManager ;
220 if (fHistEventsProcessed) delete fHistEventsProcessed ;
221 if (fCorrelation) delete fCorrelation ;
222 if (fCuts) delete fCuts;
225 //_________________________________________________________________________-
226 void AliCFTaskVertexingHF::Init()
232 if (fDebug>1) printf("AliCFTaskVertexingHF::Init()");
233 AliRDHFCuts *copyfCuts = 0x0;
235 AliFatal("No cuts defined - Exiting...");
239 switch (fDecayChannel){
241 copyfCuts = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
248 copyfCuts = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
255 copyfCuts = new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
262 copyfCuts = new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fCuts)));
269 copyfCuts = new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
276 copyfCuts = new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
279 fDauNames="K+pi+pi+pi";
283 AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exiting...");
287 const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
289 copyfCuts->SetName(nameoutput);
292 PostData(4, copyfCuts);
295 AliFatal("Failing initializing AliRDHFCuts object - Exiting...");
301 //_________________________________________________
302 void AliCFTaskVertexingHF::UserExec(Option_t *)
305 // Main loop function
308 PostData(1,fHistEventsProcessed) ;
309 PostData(2,fCFManager->GetParticleContainer()) ;
310 PostData(3,fCorrelation) ;
312 if (fFillFromGenerated){
313 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
317 Error("UserExec","NO EVENT FOUND!");
321 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
323 TClonesArray *arrayBranch=0;
325 if(!aodEvent && AODEvent() && IsStandardAOD()) {
326 // In case there is an AOD handler writing a standard AOD, use the AOD
327 // event in memory rather than the input (ESD) event.
328 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
329 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
330 // have to taken from the AOD event hold by the AliAODExtension
331 AliAODHandler* aodHandler = (AliAODHandler*)
332 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
333 if(aodHandler->GetExtensions()) {
334 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
335 AliAODEvent *aodFromExt = ext->GetAOD();
337 switch (fDecayChannel){
339 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
343 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
349 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
353 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
362 switch (fDecayChannel){
364 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
368 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
374 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
378 arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm4Prong");
386 AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
390 AliError("Could not find array of HF vertices");
396 fCFManager->SetRecEventInfo(aodEvent);
397 fCFManager->SetMCEventInfo(aodEvent);
399 //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro.
401 //loop on the MC event
403 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
405 AliError("Could not find Monte-Carlo in AOD");
410 Int_t icountReco = 0;
411 Int_t icountVertex = 0;
412 Int_t icountRefit = 0;
413 Int_t icountRecoAcc = 0;
414 Int_t icountRecoITSClusters = 0;
415 Int_t icountRecoPPR = 0;
416 Int_t icountRecoPID = 0;
419 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
421 AliError("Could not find MC Header in AOD");
425 Double_t* containerInput = new Double_t[fNvar];
426 Double_t* containerInputMC = new Double_t[fNvar];
429 AliCFVertexingHF* cfVtxHF=0x0;
430 switch (fDecayChannel){
432 cfVtxHF = new AliCFVertexingHF2Prong(mcArray, fOriginDselection);
436 cfVtxHF = new AliCFVertexingHFCascade(mcArray, fOriginDselection);
442 cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel);
446 //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection); // not there yet
453 AliError("No AliCFVertexingHF initialized");
454 delete[] containerInput;
455 delete[] containerInputMC;
459 Double_t zPrimVertex = aodVtx ->GetZ();
460 Double_t zMCVertex = mcHeader->GetVtxZ();
461 if (TMath::Abs(zMCVertex) > fCuts->GetMaxVtxZ()){
462 AliDebug(3,Form("z coordinate of MC vertex = %f, it was required to be within [-%f, +%f], skipping event", zMCVertex, fCuts->GetMaxVtxZ(), fCuts->GetMaxVtxZ()));
463 delete[] containerInput;
464 delete[] containerInputMC;
468 AliESDtrackCuts** trackCuts = new AliESDtrackCuts*[cfVtxHF->GetNProngs()];
469 if (fDecayChannel == 21){
470 // for the D*, setting the third element of the array of the track cuts to those for the soft pion
471 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs()-1; iProng++){
472 trackCuts[iProng]=fCuts->GetTrackCuts();
474 trackCuts[2] = fCuts->GetTrackCutsSoftPi();
477 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs(); iProng++){
478 trackCuts[iProng]=fCuts->GetTrackCuts();
482 //General settings: vertex, feed down and fill reco container with generated values.
483 cfVtxHF->SetRecoPrimVertex(zPrimVertex);
484 cfVtxHF->SetMCPrimaryVertex(zMCVertex);
485 cfVtxHF->SetFillFromGenerated(fFillFromGenerated);
486 cfVtxHF->SetNVar(fNvar);
487 cfVtxHF->SetFakeSelection(fFakeSelection);
488 cfVtxHF->SetRejectCandidateIfNotFromQuark(fRejectIfNoQuark);
490 // switch-off the trigger class selection (doesn't work for MC)
491 fCuts->SetTriggerClass("");
493 // MC vertex, to be used, in case, for pp
494 if (fUseMCVertex) fCuts->SetUseMCVertex();
496 if (fCentralitySelection){ // keep only the requested centrality
497 if(fCuts->IsEventSelectedInCentrality(aodEvent)!=0) {
498 delete[] containerInput;
499 delete[] containerInputMC;
503 } else { // keep all centralities
504 fCuts->SetMinCentrality(0.);
505 fCuts->SetMaxCentrality(100.);
509 Float_t centValue = fCuts->GetCentrality(aodEvent);
510 cfVtxHF->SetCentralityValue(centValue);
512 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
513 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
515 AliError("Failed casting particle from MC array!, Skipping particle");
518 // check the MC-level cuts, must be the desidered particle
519 if (!fCFManager->CheckParticleCuts(0, mcPart)) {
520 continue; // 0 stands for MC level
522 cfVtxHF->SetMCCandidateParam(iPart);
525 cquarks += cfVtxHF->MCcquarkCounting(mcPart);
527 if (!(cfVtxHF->SetLabelArray())){
528 AliDebug(2,Form("Impossible to set the label array (decaychannel = %d)",fDecayChannel));
532 //check the candiate family at MC level
533 if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) {
534 AliDebug(2,Form("Check on the family wrong!!! (decaychannel = %d)",fDecayChannel));
538 AliDebug(2,Form("Check on the family OK!!! (decaychannel = %d)",fDecayChannel));
541 //Fill the MC container
542 Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
543 if (mcContainerFilled) {
544 if (fUseWeight)fWeight = GetWeight(containerInputMC[0]);
545 if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) continue;
546 //MC Limited Acceptance
547 if (TMath::Abs(containerInputMC[1]) < 0.5) {
548 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight);
549 AliDebug(3,"MC Lim Acc container filled\n");
553 fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight);
555 AliDebug(3,"MC cointainer filled \n");
558 // check the MC-Acceptance level cuts
559 // since standard CF functions are not applicable, using Kine Cuts on daughters
560 Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep();
562 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight);
563 AliDebug(3,"MC acceptance cut passed\n");
567 if (fCuts->IsEventSelected(aodEvent)){
568 // filling the container if the vertex is ok
569 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex, fWeight) ;
570 AliDebug(3,"Vertex cut passed and container filled\n");
573 //mc Refit requirement
574 Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, &trackCuts[0]);
576 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight);
577 AliDebug(3,"MC Refit cut passed and container filled\n");
583 AliDebug(3,"MC Refit cut not passed\n");
588 AliDebug (3, "MC vertex step not passed\n");
593 AliDebug (3, "MC in acceptance step not passed\n");
598 AliDebug (3, "MC container not filled\n");
602 if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks));
603 AliDebug(2,Form("Found %i MC particles that are %s!!",icountMC,fPartName.Data()));
604 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Acc cuts!!",icountAcc,fPartName.Data()));
605 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Vertex cuts!!",icountVertex,fPartName.Data()));
606 AliDebug(2,Form("Found %i MC particles that are %s and satisfy Refit cuts!!",icountRefit,fPartName.Data()));
608 // Now go to rec level
609 fCountMC += icountMC;
610 fCountAcc += icountAcc;
611 fCountVertex+= icountVertex;
612 fCountRefit+= icountRefit;
614 AliDebug(2,Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel));
616 for(Int_t iCandid = 0; iCandid<arrayBranch->GetEntriesFast();iCandid++){
617 AliAODRecoDecayHF* charmCandidate=0x0;
618 switch (fDecayChannel){
620 charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid);
624 charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid);
630 charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid);
634 charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid);
641 Bool_t unsetvtx=kFALSE;
642 if(!charmCandidate->GetOwnPrimaryVtx()) {
643 charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
647 Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate);
648 if (!signAssociation){
649 charmCandidate = 0x0;
653 Int_t isPartOrAntipart = cfVtxHF->CheckReflexion(fSign);
654 if (isPartOrAntipart == 0){
655 AliDebug(2, Form("The candidate pdg code doesn't match the requirement set in the task (fSign = %d)",fSign));
660 Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
663 // weight according to pt
664 if (fUseWeight)fWeight = GetWeight(containerInput[0]);
666 if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue;
669 Bool_t recoStep = cfVtxHF->RecoStep();
670 Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
672 if (recoStep && recoContFilled && vtxCheck){
673 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
675 AliDebug(3,"Reco step passed and container filled\n");
677 //Reco in the acceptance -- take care of UNFOLDING!!!!
678 Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(&trackCuts[0]);
679 if (recoAcceptanceStep) {
680 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
682 AliDebug(3,"Reco acceptance cut passed and container filled\n");
685 Double_t fill[4]; //fill response matrix
686 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
687 if (bUnfolding) fCorrelation->Fill(fill);
690 //Number of ITS cluster requirements
691 Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
692 if (recoITSnCluster){
693 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
694 icountRecoITSClusters++;
695 AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
697 Bool_t iscutsusingpid = fCuts->GetIsUsePID();
698 Int_t recoAnalysisCuts = -1, recoPidSelection = -1;
699 fCuts->SetUsePID(kFALSE);
700 recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
702 if (fDecayChannel==33){ // Ds case, where more possibilities are considered
703 Bool_t keepDs=ProcessDs(recoAnalysisCuts);
704 if(keepDs) recoAnalysisCuts=3;
708 fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object
709 if (recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart){
710 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight);
712 AliDebug(3,"Reco Analysis cuts passed and container filled \n");
714 //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID);
715 //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){
716 recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
718 if (fDecayChannel==33){ // Ds case, where more possibilities are considered
719 Bool_t keepDs=ProcessDs(recoPidSelection);
720 if(keepDs) recoPidSelection=3;
723 if (recoPidSelection == 3 || recoPidSelection == isPartOrAntipart){
724 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight);
726 AliDebug(3,"Reco PID cuts passed and container filled \n");
728 Double_t fill[4]; //fill response matrix
729 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
730 if (bUnfolding) fCorrelation->Fill(fill);
734 AliDebug(3, "Analysis Cuts step not passed \n");
739 AliDebug(3, "PID selection not passed \n");
744 AliDebug(3, "Number of ITS cluster step not passed\n");
749 AliDebug(3, "Reco acceptance step not passed\n");
754 AliDebug(3, "Reco step not passed\n");
759 if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
760 } // end loop on candidate
762 fCountReco+= icountReco;
763 fCountRecoAcc+= icountRecoAcc;
764 fCountRecoITSClusters+= icountRecoITSClusters;
765 fCountRecoPPR+= icountRecoPPR;
766 fCountRecoPID+= icountRecoPID;
768 fHistEventsProcessed->Fill(0);
770 delete[] containerInput;
771 delete[] containerInputMC;
774 // for (Int_t i=0; i<cfVtxHF->GetNProngs(); i++){
775 // delete [] trackCuts[i];
781 //___________________________________________________________________________
782 void AliCFTaskVertexingHF::Terminate(Option_t*)
784 // The Terminate() function is the last function to be called during
785 // a query. It always runs on the client, it can be used to present
786 // the results graphically or save the results to file.
788 AliAnalysisTaskSE::Terminate();
790 AliInfo(Form("Found %i MC particles that are %s in MC, in %d events",fCountMC,fPartName.Data(),fEvents));
791 AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, in %d events",fCountAcc,fPartName.Data(),fEvents));
792 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));
793 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));
794 AliInfo(Form("Found %i reco %s that are decaying in %s, in %d events",fCountReco,fPartName.Data(),fDauNames.Data(),fEvents));
795 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));
796 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));
797 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));
798 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));
800 // draw some example plots....
802 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
804 printf("CONTAINER NOT FOUND\n");
807 // projecting the containers to obtain histograms
808 // first argument = variable, second argument = step
811 for(Int_t iC=0;iC<12; iC++){
813 h[0][iC] = cont->ShowProjection(iC,0);
814 // MC-Acceptance level
815 h[1][iC] = cont->ShowProjection(iC,1);
817 h[2][iC] = cont->ShowProjection(iC,4);
820 if(fDecayChannel==31){
821 titles[0]="pT_Dplus (GeV/c)";
822 titles[1]="rapidity";
823 titles[2]="phi (rad)";
824 titles[3]="cT (#mum)";
825 titles[4]="cosPointingAngle";
826 titles[5]="pT_1 (GeV/c)";
827 titles[6]="pT_2 (GeV/c)";
828 titles[7]="pT_3 (GeV/c)";
829 titles[8]="d0_1 (#mum)";
830 titles[9]="d0_2 (#mum)";
831 titles[10]="d0_3 (#mum)";
832 titles[11]="zVertex (cm)";
834 titles[0]="pT_D0 (GeV/c)";
835 titles[1]="rapidity";
836 titles[2]="cosThetaStar";
837 titles[3]="pT_pi (GeV/c)";
838 titles[4]="pT_K (Gev/c)";
839 titles[5]="cT (#mum)";
840 titles[6]="dca (#mum)";
841 titles[7]="d0_pi (#mum)";
842 titles[8]="d0_K (#mum)";
843 titles[9]="d0xd0 (#mum^2)";
844 titles[10]="cosPointingAngle";
845 titles[11]="phi (rad)";
847 Int_t markers[12]={20,24,21,25,27,28,
849 Int_t colors[3]={2,8,4};
850 for(Int_t iC=0;iC<12; iC++){
851 for(Int_t iStep=0;iStep<3;iStep++){
852 h[iStep][iC]->SetTitle(titles[iC].Data());
853 h[iStep][iC]->GetXaxis()->SetTitle(titles[iC].Data());
854 Double_t maxh=h[iStep][iC]->GetMaximum();
855 h[iStep][iC]->GetYaxis()->SetRangeUser(0,maxh*1.2);
856 h[iStep][iC]->SetMarkerStyle(markers[iC]);
857 h[iStep][iC]->SetMarkerColor(colors[iStep]);
863 gStyle->SetCanvasColor(0);
864 gStyle->SetFrameFillColor(0);
865 gStyle->SetTitleFillColor(0);
866 gStyle->SetStatColor(0);
868 // drawing in 2 separate canvas for a matter of clearity
869 TCanvas * c1 =new TCanvas("c1New","Vars 0,1,2",1100,1600);
872 for(Int_t iVar=0; iVar<3; iVar++){
874 h[0][iVar]->Draw("p");
876 h[1][iVar]->Draw("p");
878 h[2][iVar]->Draw("p");
881 TCanvas * c2 =new TCanvas("c2New","Vars 3,4,5",1100,1600);
884 for(Int_t iVar=3; iVar<6; iVar++){
886 h[0][iVar]->Draw("p");
888 h[1][iVar]->Draw("p");
890 h[2][iVar]->Draw("p");
894 TCanvas * c3 =new TCanvas("c3New","Vars 6,7,8",1100,1600);
897 for(Int_t iVar=6; iVar<9; iVar++){
899 h[0][iVar]->Draw("p");
901 h[1][iVar]->Draw("p");
903 h[2][iVar]->Draw("p");
907 TCanvas * c4 =new TCanvas("c4New","Vars 9,10,11",1100,1600);
910 for(Int_t iVar=9; iVar<11; iVar++){
912 h[0][iVar]->Draw("p");
914 h[1][iVar]->Draw("p");
916 h[2][iVar]->Draw("p");
920 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
922 TH2D* corr1 =hcorr->Projection(0,2);
923 TH2D* corr2 = hcorr->Projection(1,3);
925 TCanvas * c7 =new TCanvas("c7New","",800,400);
933 TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
937 for(Int_t iC=0;iC<12; iC++){
938 for(Int_t iStep=0;iStep<3;iStep++){
939 h[iStep][iC]->Write(Form("Step%d_%s",iStep,titles[iC].Data()));
942 file_projection->Close();
947 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
948 c2->SaveAs("Plots/pTpi_pTK_cT.eps");
949 c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
950 c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
952 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
953 c2->SaveAs("Plots/pTpi_pTK_cT.gif");
954 c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
955 c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
959 //___________________________________________________________________________
960 void AliCFTaskVertexingHF::UserCreateOutputObjects()
962 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
963 //TO BE SET BEFORE THE EXECUTION OF THE TASK
965 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
969 fHistEventsProcessed = new TH1I("CFHFchist0","",1,0,1) ;
971 PostData(1,fHistEventsProcessed) ;
972 PostData(2,fCFManager->GetParticleContainer()) ;
973 PostData(3,fCorrelation) ;
978 //_________________________________________________________________________
979 Double_t AliCFTaskVertexingHF::GetWeight(Float_t pt)
982 // calculating the weight to fill the container
986 // p0 = 1.63297e-01 --> 0.322643
992 // p0 = 1.85906e-01 --> 0.36609
997 Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
998 Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
1000 Double_t dndpt_func1 = dNdptFit(pt,func1);
1001 Double_t dndpt_func2 = dNdptFit(pt,func2);
1002 AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
1003 return dndpt_func1/dndpt_func2;
1006 //__________________________________________________________________________________________________
1007 Double_t AliCFTaskVertexingHF::dNdptFit(Float_t pt, Double_t* par)
1010 // calculating dNdpt
1013 Double_t denom = TMath::Power((pt/par[1]), par[3] );
1014 Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
1020 //__________________________________________________________________________________________________
1021 Bool_t AliCFTaskVertexingHF::ProcessDs(Int_t recoAnalysisCuts) const{
1022 // processes output of Ds is selected
1024 if(recoAnalysisCuts > 0){
1025 Int_t isKKpi=recoAnalysisCuts&1;
1026 Int_t ispiKK=recoAnalysisCuts&2;
1027 Int_t isPhiKKpi=recoAnalysisCuts&4;
1028 Int_t isPhipiKK=recoAnalysisCuts&8;
1029 Int_t isK0starKKpi=recoAnalysisCuts&16;
1030 Int_t isK0starpiKK=recoAnalysisCuts&32;
1032 if(isKKpi && isPhiKKpi) keep=kTRUE;
1033 if(ispiKK && isPhipiKK) keep=kTRUE;
1035 else if(fDsOption==2){
1036 if(isKKpi && isK0starKKpi) keep=kTRUE;
1037 if(ispiKK && isK0starpiKK) keep=kTRUE;
1039 else if(fDsOption==3)keep=kTRUE;