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 **************************************************************************/
18 //-----------------------------------------------------------------------
19 // Class for HF corrections as a function of many variables
20 // 6 Steps introduced: MC, MC Acc, Reco, Reco Acc, Reco Acc + ITS Cl,
21 // Reco Acc + ITS Cl + PPR cuts
22 // 13 variables used: pt, y, cosThetaStar, ptPi, ptK, ct,
23 // dca, d0Pi, d0K, d0Pixd0K, cosPointingAngle, phi, z
25 //-----------------------------------------------------------------------
26 // Author : C. Zampolli, CERN
27 //-----------------------------------------------------------------------
28 //-----------------------------------------------------------------------
29 // Base class for HF Unfolding (pt and eta)
30 // correlation matrix filled at Acceptance and PPR level
31 // Author: A.Grelli , Utrecht - agrelli@uu.nl
32 //-----------------------------------------------------------------------
34 #include <TParticle.h>
35 #include <TDatabasePDG.h>
40 #include "AliCFHeavyFlavourTaskMultiVarMultiStep.h"
42 #include "AliMCEvent.h"
43 #include "AliCFManager.h"
44 #include "AliCFContainer.h"
46 #include "AliAnalysisManager.h"
47 #include "AliAODHandler.h"
48 #include "AliAODEvent.h"
49 #include "AliAODRecoDecay.h"
50 #include "AliAODRecoDecayHF.h"
51 #include "AliAODRecoDecayHF2Prong.h"
52 #include "AliAODMCParticle.h"
53 #include "AliAODMCHeader.h"
54 #include "AliESDtrack.h"
55 #include "AliRDHFCutsD0toKpi.h"
57 #include "THnSparse.h"
59 #include "AliAnalysisDataSlot.h"
60 #include "AliAnalysisDataContainer.h"
62 //__________________________________________________________________________
63 AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep() :
67 fHistEventsProcessed(0x0),
75 fCountRecoITSClusters(0),
79 fFillFromGenerated(kFALSE),
81 fAcceptanceUnf(kTRUE),
83 fKeepD0fromBOnly(kFALSE),
93 //___________________________________________________________________________
94 AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(const Char_t* name, AliRDHFCutsD0toKpi* cuts) :
95 AliAnalysisTaskSE(name),
98 fHistEventsProcessed(0x0),
106 fCountRecoITSClusters(0),
110 fFillFromGenerated(kFALSE),
112 fAcceptanceUnf(kTRUE),
113 fKeepD0fromB(kFALSE),
114 fKeepD0fromBOnly(kFALSE),
121 // Constructor. Initialization of Inputs and Outputs
123 Info("AliCFHeavyFlavourTaskMultiVarMultiStep","Calling Constructor");
125 DefineInput(0) and DefineOutput(0)
126 are taken care of by AliAnalysisTaskSE constructor
128 DefineOutput(1,TH1I::Class());
129 DefineOutput(2,AliCFContainer::Class());
130 DefineOutput(3,THnSparseD::Class());
131 DefineOutput(4,AliRDHFCutsD0toKpi::Class());
136 //___________________________________________________________________________
137 AliCFHeavyFlavourTaskMultiVarMultiStep& AliCFHeavyFlavourTaskMultiVarMultiStep::operator=(const AliCFHeavyFlavourTaskMultiVarMultiStep& c)
140 // Assignment operator
143 AliAnalysisTaskSE::operator=(c) ;
145 fCFManager = c.fCFManager;
146 fHistEventsProcessed = c.fHistEventsProcessed;
152 //___________________________________________________________________________
153 AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(const AliCFHeavyFlavourTaskMultiVarMultiStep& c) :
154 AliAnalysisTaskSE(c),
156 fCFManager(c.fCFManager),
157 fHistEventsProcessed(c.fHistEventsProcessed),
158 fCorrelation(c.fCorrelation),
159 fCountMC(c.fCountMC),
160 fCountAcc(c.fCountAcc),
161 fCountVertex(c.fCountVertex),
162 fCountRefit(c.fCountRefit),
163 fCountReco(c.fCountReco),
164 fCountRecoAcc(c.fCountRecoAcc),
165 fCountRecoITSClusters(c.fCountRecoITSClusters),
166 fCountRecoPPR(c.fCountRecoPPR),
167 fCountRecoPID(c.fCountRecoPID),
169 fFillFromGenerated(c.fFillFromGenerated),
170 fMinITSClusters(c.fMinITSClusters),
171 fAcceptanceUnf(c.fAcceptanceUnf),
172 fKeepD0fromB(c.fKeepD0fromB),
173 fKeepD0fromBOnly(c.fKeepD0fromBOnly),
175 fUseWeight(c.fUseWeight),
184 //___________________________________________________________________________
185 AliCFHeavyFlavourTaskMultiVarMultiStep::~AliCFHeavyFlavourTaskMultiVarMultiStep() {
189 if (fCFManager) delete fCFManager ;
190 if (fHistEventsProcessed) delete fHistEventsProcessed ;
191 if (fCorrelation) delete fCorrelation ;
192 // if (fCuts) delete fCuts ;
195 //________________________________________________________________________
196 void AliCFHeavyFlavourTaskMultiVarMultiStep::Init(){
202 if(fDebug > 1) printf("AliCFHeavyFlavourTaskMultiVarMultiStep::Init() \n");
204 fMinITSClusters = fCuts->GetTrackCuts()->GetMinNClustersITS();
205 AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*fCuts);
206 const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
207 copyfCuts->SetName(nameoutput);
209 PostData(4,copyfCuts);
214 //_________________________________________________
215 void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
218 // Main loop function
221 PostData(1,fHistEventsProcessed) ;
222 PostData(2,fCFManager->GetParticleContainer()) ;
223 PostData(3,fCorrelation) ;
225 AliESDtrackCuts* trackCuts = fCuts->GetTrackCuts(); // track cuts
227 if (fFillFromGenerated){
228 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
232 Error("UserExec","NO EVENT FOUND!");
236 // check that the fKeepD0fromB flag is set to true when the fKeepD0fromBOnly flag is true
237 if(fKeepD0fromBOnly) {
239 if(fEvents<2) AliInfo(Form("Both fKeepD0fromB and fKeepD0fromBOnly flags are true, looking _ONLY_ at D0 FROM B"));
242 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
244 TClonesArray *arrayD0toKpi=0;
246 if(!aodEvent && AODEvent() && IsStandardAOD()) {
247 // In case there is an AOD handler writing a standard AOD, use the AOD
248 // event in memory rather than the input (ESD) event.
249 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
250 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
251 // have to taken from the AOD event hold by the AliAODExtension
252 AliAODHandler* aodHandler = (AliAODHandler*)
253 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
254 if(aodHandler->GetExtensions()) {
255 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
256 AliAODEvent *aodFromExt = ext->GetAOD();
257 arrayD0toKpi=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
260 arrayD0toKpi=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
265 AliError("Could not find array of HF vertices");
269 // fix for temporary bug in ESDfilter
270 // the AODs with null vertex pointer didn't pass the PhysSel
271 if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
275 fCFManager->SetRecEventInfo(aodEvent);
276 fCFManager->SetMCEventInfo(aodEvent);
278 // MC-event selection
279 Double_t containerInput[13] ;
280 Double_t containerInputMC[13] ;
282 //loop on the MC event
284 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
286 AliError("Could not find Monte-Carlo in AOD");
291 Int_t icountReco = 0;
292 Int_t icountVertex = 0;
293 Int_t icountRefit = 0;
294 Int_t icountRecoAcc = 0;
295 Int_t icountRecoITSClusters = 0;
296 Int_t icountRecoPPR = 0;
297 Int_t icountRecoPID = 0;
299 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
301 AliError("Could not find MC Header in AOD");
307 // AOD primary vertex
308 AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
310 AliError("There is no primary vertex !");
313 Double_t zPrimVertex = vtx1->GetZ();
314 Double_t zMCVertex = mcHeader->GetVtxZ();
315 Bool_t vtxFlag = kTRUE;
316 TString title=vtx1->GetTitle();
317 if(!title.Contains("VertexerTracks")) vtxFlag=kFALSE;
319 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
320 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
322 AliWarning("Particle not found in tree, skipping");
325 if (mcPart->GetPdgCode() == 4) cquarks++;
326 if (mcPart->GetPdgCode() == -4) cquarks++;
328 // check the MC-level cuts
330 if (!fCFManager->CheckParticleCuts(0, mcPart)) continue; // 0 stands for MC level
331 Int_t pdgGranma = CheckOrigin(mcPart, mcArray);
332 Int_t abspdgGranma = TMath::Abs(pdgGranma);
333 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
334 AliDebug(2,Form("Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
335 if (!fKeepD0fromB) continue; // skipping particles that don't come from c quark
337 else { if(fKeepD0fromBOnly) continue; } // skipping particles that don't come from b quark
339 // if (TMath::Abs(pdgGranma)!=4) {
341 // fill the container for Gen-level selection
342 Double_t vectorMC[7] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.};
344 if (GetGeneratedValuesFromMCParticle(mcPart, mcArray, vectorMC)){
345 containerInputMC[0] = vectorMC[0];
346 containerInputMC[1] = vectorMC[1] ;
347 containerInputMC[2] = vectorMC[2] ;
348 containerInputMC[3] = vectorMC[3] ;
349 containerInputMC[4] = vectorMC[4] ;
350 containerInputMC[5] = vectorMC[5] ; // in micron
351 containerInputMC[6] = 0.; // dummy value, meaningless in MC, in micron
352 containerInputMC[7] = 0.; // dummy value, meaningless in MC, in micron
353 containerInputMC[8] = 0.; // dummy value, meaningless in MC, in micron
354 containerInputMC[9] = -100000.; // dummy value, meaningless in MC, in micron^2
355 containerInputMC[10] = 1.01; // dummy value, meaningless in MC
356 containerInputMC[11] = vectorMC[6]; // dummy value, meaningless in MC
357 containerInputMC[12] = zMCVertex; // z of reconstructed of primary vertex
358 if (fUseWeight) fWeight = GetWeight(vectorMC[0]); // setting the weight according to the function defined in AliCFHeavyFlavourTaskMultiVarMultiStep::GetWeight(Float_t pt)
359 AliDebug(3,Form("weight = %f",fWeight));
360 if (!fCuts->IsInFiducialAcceptance(vectorMC[0],vectorMC[1])) continue;
361 if (TMath::Abs(vectorMC[1]) < 0.5) {
362 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc,fWeight);
364 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGenerated,fWeight);
367 // check the MC-Acceptance level cuts
368 // since standard CF functions are not applicable, using Kine Cuts on daughters
370 Int_t daughter0 = mcPart->GetDaughter(0);
371 Int_t daughter1 = mcPart->GetDaughter(1);
372 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
373 if (daughter0 == 0 || daughter1 == 0) {
374 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!! But it should have, this check was already done...");
376 if (TMath::Abs(daughter1 - daughter0) != 1) {
377 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, but it should be, this check was already done...");
379 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
380 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
381 if (!mcPartDaughter0 || !mcPartDaughter1) {
382 AliWarning("At least one Daughter Particle not found in tree, but it should be, this check was already done...");
385 Double_t eta0 = mcPartDaughter0->Eta();
386 Double_t eta1 = mcPartDaughter1->Eta();
387 Double_t y0 = mcPartDaughter0->Y();
388 Double_t y1 = mcPartDaughter1->Y();
389 Double_t pt0 = mcPartDaughter0->Pt();
390 Double_t pt1 = mcPartDaughter1->Pt();
391 AliDebug(2, Form("Daughter 0: eta = %f, y = %f, pt = %f", eta0, y0, pt0));
392 AliDebug(2, Form("Daughter 1: eta = %f, y = %f, pt = %f", eta1, y1, pt1));
393 Bool_t daught0inAcceptance = (TMath::Abs(eta0) < 0.9 && pt0 > 0.1);
394 Bool_t daught1inAcceptance = (TMath::Abs(eta1) < 0.9 && pt1 > 0.1);
395 if (daught0inAcceptance && daught1inAcceptance) {
396 // checking whether the cuts implemented in the CF are equivalent - simply a cross-check
397 AliDebug(2, "Daughter particles in acceptance");
398 if (!fCFManager->CheckParticleCuts(1, mcPartDaughter0)) {
399 AliDebug(2,"Inconsistency with CF cut for daughter 0!");
401 if (!fCFManager->CheckParticleCuts(1, mcPartDaughter1)) {
402 AliDebug(2,"Inconsistency with CF cut for daughter 1!");
404 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance,fWeight);
407 // check on the vertex
408 if (fCuts->IsEventSelected(aodEvent)){
409 AliDebug(2,"Vertex cut passed\n");
410 // filling the container if the vertex is ok
411 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex,fWeight) ;
413 // check on the kTPCrefit and kITSrefit conditions of the daughters
414 Bool_t refitFlag = kTRUE;
415 if (trackCuts->GetRequireTPCRefit() || trackCuts->GetRequireITSRefit()){
416 Int_t foundDaughters = 0;
417 for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
418 AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
419 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
420 if ((track->GetLabel() == daughter0) || (track->GetLabel() == daughter1)) {
422 if (trackCuts->GetRequireTPCRefit()) {
423 if(!(track->GetStatus()&AliESDtrack::kTPCrefit)){
428 if (trackCuts->GetRequireITSRefit()) {
429 if(!(track->GetStatus()&AliESDtrack::kITSrefit)){
435 if (foundDaughters == 2){ // both daughters have been checked
439 if (foundDaughters != 2) refitFlag = kFALSE;
442 AliDebug(3,"Refit cut passed\n");
443 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit,fWeight);
447 AliDebug(3,"Refit cut not passed\n");
451 AliDebug(3,"Vertex cut not passed\n");
455 AliDebug(3,"Acceptance cut not passed\n");
459 AliDebug(3,"Problems in filling the container");
464 if (cquarks<2) AliDebug(2,Form("Event found with %d c-quarks", cquarks));
466 AliDebug(2,Form("Found %i MC particles that are D0!!",icountMC));
467 AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Acc cuts!!",icountAcc));
468 AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Vertex cuts!!",icountVertex));
469 AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Refit cuts!!",icountRefit));
471 // Now go to rec level
472 fCountMC += icountMC;
473 fCountAcc += icountAcc;
475 AliDebug(2, Form("Found %d vertices",arrayD0toKpi->GetEntriesFast()));
477 Int_t pdgDgD0toKpi[2]={321,211};
478 Int_t isD0D0bar=1;// 1 for D0, 2 for D0bar, to be used for the PPR and PID selection steps
480 for (Int_t iD0toKpi = 0; iD0toKpi<arrayD0toKpi->GetEntriesFast(); iD0toKpi++) {
482 AliAODRecoDecayHF2Prong* d0tokpi = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->At(iD0toKpi);
483 if(!d0tokpi) continue;
484 Bool_t unsetvtx=kFALSE;
485 if(!d0tokpi->GetOwnPrimaryVtx()) {
486 d0tokpi->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
490 // find associated MC particle
491 Int_t mcLabel = d0tokpi->MatchToMC(421,mcArray,2,pdgDgD0toKpi) ;
494 AliDebug(2,"No MC particle found");
498 AliAODMCParticle* mcVtxHF = (AliAODMCParticle*)mcArray->At(mcLabel);
500 AliWarning("Could not find associated MC in AOD MC tree");
504 if (mcVtxHF->GetPdgCode() == 421){ // particle is D0
505 if (fSign == 1){ // I ask for D0bar only
506 AliDebug(2,"particle is D0, I ask for D0bar only");
513 else if (mcVtxHF->GetPdgCode()== -421){ // particle is D0bar
514 if (fSign == 0){ // I ask for D0 only
515 AliDebug(2,"particle is D0bar, I ask for D0 only");
524 // check whether the daughters have kTPCrefit and kITSrefit set
525 AliAODTrack *track0 = (AliAODTrack*)d0tokpi->GetDaughter(0);
526 AliAODTrack *track1 = (AliAODTrack*)d0tokpi->GetDaughter(1);
527 if( !track0 || !track1 ) {
528 AliWarning("Could not find associated MC daughter tracks in AOD MC tree");
531 if ((trackCuts->GetRequireTPCRefit() && (!(track0->GetStatus()&AliESDtrack::kTPCrefit) || !(track1->GetStatus()&AliESDtrack::kTPCrefit))) ||
532 (trackCuts->GetRequireITSRefit() && (!(track0->GetStatus()&AliESDtrack::kITSrefit) || !(track1->GetStatus()&AliESDtrack::kITSrefit)))){
533 // skipping if at least one daughter does not have kTPCrefit or kITSrefit, if they were required
537 // check on the vertex -- could be moved outside the loop on the reconstructed D0...
538 if(!fCuts->IsEventSelected(aodEvent)) {
539 // skipping cause vertex is not ok
543 const Double_t d0tokpiCuts[9] = {0.3,999999.,1.1,0.,0.,999999.,999999.,999999.,0.};
545 if (!(d0tokpi->SelectD0(&d0tokpiCuts[0],okD0,okD0bar))){
546 // skipping candidate
550 // check if associated MC v0 passes the cuts
551 if (!fCFManager->CheckParticleCuts(0 ,mcVtxHF)) { // 0 stands for MC
552 AliDebug(2, "Skipping the particles due to cuts");
555 Int_t pdgGranma = CheckOrigin(mcVtxHF, mcArray);
556 Int_t abspdgGranma = TMath::Abs(pdgGranma);
557 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
558 AliDebug(2,Form("At Reco level, from MC info: Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
559 if (!fKeepD0fromB) continue; // skipping particles that don't come from c quark
561 else { if(fKeepD0fromBOnly) continue; } // skipping particles that don't come from b quark
563 // fill the container...
564 Double_t pt = d0tokpi->Pt();
565 Double_t rapidity = d0tokpi->YD0();
567 Double_t cosThetaStar = 9999.;
570 Double_t dca = d0tokpi->GetDCA();
573 Double_t d0xd0 = d0tokpi->Prodd0d0();
574 Double_t cosPointingAngle = d0tokpi->CosPointingAngle();
575 Double_t phi = d0tokpi->Phi();
576 Int_t pdgCode = mcVtxHF->GetPdgCode();
578 cosThetaStar = d0tokpi->CosThetaStarD0();
579 pTpi = d0tokpi->PtProng(0);
580 pTK = d0tokpi->PtProng(1);
581 d0pi = d0tokpi->Getd0Prong(0);
582 d0K = d0tokpi->Getd0Prong(1);
583 invMass=d0tokpi->InvMassD0();
586 cosThetaStar = d0tokpi->CosThetaStarD0bar();
587 pTpi = d0tokpi->PtProng(1);
588 pTK = d0tokpi->PtProng(0);
589 d0pi = d0tokpi->Getd0Prong(1);
590 d0K = d0tokpi->Getd0Prong(0);
591 invMass=d0tokpi->InvMassD0bar();
594 Double_t cT = d0tokpi->CtD0();
596 if (!fFillFromGenerated){
597 // ...either with reconstructed values....
598 containerInput[0] = pt;
599 containerInput[1] = rapidity;
600 containerInput[2] = cosThetaStar;
601 containerInput[3] = pTpi;
602 containerInput[4] = pTK;
603 containerInput[5] = cT*1.E4; // in micron
604 containerInput[6] = dca*1.E4; // in micron
605 containerInput[7] = d0pi*1.E4; // in micron
606 containerInput[8] = d0K*1.E4; // in micron
607 containerInput[9] = d0xd0*1.E8; // in micron^2
608 containerInput[10] = cosPointingAngle; // in micron
609 containerInput[11] = phi;
610 containerInput[12] = zPrimVertex; // z of reconstructed of primary vertex
613 // ... or with generated values
614 Double_t vectorMC[7] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.};
615 if (GetGeneratedValuesFromMCParticle(mcVtxHF, mcArray, vectorMC)){
616 containerInput[0] = vectorMC[0];
617 containerInput[1] = vectorMC[1] ;
618 containerInput[2] = vectorMC[2] ;
619 containerInput[3] = vectorMC[3] ;
620 containerInput[4] = vectorMC[4] ;
621 containerInput[5] = vectorMC[5] ; // in micron
622 containerInput[6] = 0.; // dummy value, meaningless in MC, in micron
623 containerInput[7] = 0.; // dummy value, meaningless in MC, in micron
624 containerInput[8] = 0.; // dummy value, meaningless in MC, in micron
625 containerInput[9] = 100000.; // dummy value, meaningless in MC, in micron^2
626 containerInput[10] = 1.01; // dummy value, meaningless in MC
627 containerInput[11] = vectorMC[6];
628 containerInput[12] = zMCVertex; // z of reconstructed of primary vertex
631 AliDebug(3,"Problems in filling the container");
636 if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue; // fiducial region
638 AliDebug(2, Form("Filling the container with pt = %f, rapidity = %f, cosThetaStar = %f, pTpi = %f, pTK = %f, cT = %f", containerInput[0], containerInput[1], containerInput[2], containerInput[3], containerInput[4], containerInput[5]));
640 AliDebug(2,Form("%d: filling RECO step\n",iD0toKpi));
641 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed,fWeight) ;
644 Float_t etaCutMin, etaCutMax, ptCutMin, ptCutMax;
645 trackCuts->GetEtaRange(etaCutMin,etaCutMax);
646 trackCuts->GetPtRange(ptCutMin,ptCutMax);
647 Bool_t acceptanceProng0 = (d0tokpi->EtaProng(0)>etaCutMin && d0tokpi->EtaProng(0)<etaCutMax && d0tokpi->PtProng(0) > ptCutMin && d0tokpi->PtProng(0) < ptCutMax);
648 Bool_t acceptanceProng1 = (d0tokpi->EtaProng(1)>etaCutMin && d0tokpi->EtaProng(1)<etaCutMax && d0tokpi->PtProng(1) > ptCutMin && d0tokpi->PtProng(1) < ptCutMax);
649 if (acceptanceProng0 && acceptanceProng1) {
650 AliDebug(2,"D0 reco daughters in acceptance");
651 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance,fWeight) ;
656 Double_t fill[4]; //fill response matrix
658 // dimensions 0&1 : pt,eta (Rec)
663 // dimensions 2&3 : pt,eta (MC)
665 fill[2] = mcVtxHF->Pt();
666 fill[3] = mcVtxHF->Y();
668 fCorrelation->Fill(fill);
672 // cut on the min n. of clusters in ITS
673 if (fCuts->IsSelected(d0tokpi,AliRDHFCuts::kTracks)){
674 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
675 icountRecoITSClusters++;
676 AliDebug(2,Form("pT = %f, dca = %f, cosThetaStar = %f, pTpi = %f, pTK = %f, d0pi = %f, d0K = %f, d0xd0 = %f, cosPointingAngle = %f", pt, dca, cosThetaStar,pTpi, pTK, d0pi*1E4, d0K*1E4, d0xd0*1E8, cosPointingAngle));
678 // setting the use of the PID cut when applying the selection to FALSE - whatever it was. Keeping track of the original value
679 Bool_t iscutusingpid=fCuts->GetIsUsePID();
680 Int_t isselcuts=-1,isselpid=-1;
681 fCuts->SetUsePID(kFALSE);
682 //Bool_t origFlag = fCuts->GetIsPrimaryWithoutDaughters();
683 //fCuts->SetRemoveDaughtersFromPrim(kFALSE);
684 isselcuts = fCuts->IsSelected(d0tokpi,AliRDHFCuts::kCandidate,aodEvent);
685 //fCuts->SetRemoveDaughtersFromPrim(origFlag);
686 fCuts->SetUsePID(iscutusingpid); // restoring usage of the PID from the cuts object
687 if (isselcuts == 3 || isselcuts == isD0D0bar){
688 AliDebug(2,"Particle passed PPR cuts (actually cuts for D0 analysis!)");
689 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPPR,fWeight) ;
692 if(!fAcceptanceUnf){ // unfolding
694 Double_t fill[4]; //fill response matrix
696 // dimensions 0&1 : pt,eta (Rec)
701 // dimensions 2&3 : pt,eta (MC)
703 fill[2] = mcVtxHF->Pt();
704 fill[3] = mcVtxHF->Y();
706 fCorrelation->Fill(fill);
710 isselpid = fCuts->IsSelected(d0tokpi,AliRDHFCuts::kPID);
711 if((fCuts->CombineSelectionLevels(3,isselcuts,isselpid)==isD0D0bar)||(fCuts->CombineSelectionLevels(3,isselcuts,isselpid)==3)){
712 AliDebug(2,"Particle passed PID cuts");
713 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPID,fWeight) ;
720 if(unsetvtx) d0tokpi->UnsetOwnPrimaryVtx();
721 } // end loop on D0->Kpi
723 AliDebug(2, Form("Found %i Reco particles that are D0!!",icountReco));
725 fCountReco+= icountReco;
726 fCountVertex+= icountVertex;
727 fCountRefit+= icountRefit;
728 fCountRecoAcc+= icountRecoAcc;
729 fCountRecoITSClusters+= icountRecoITSClusters;
730 fCountRecoPPR+= icountRecoPPR;
731 fCountRecoPID+= icountRecoPID;
733 fHistEventsProcessed->Fill(0);
734 /* PostData(0) is taken care of by AliAnalysisTaskSE */
735 //PostData(1,fHistEventsProcessed) ;
736 //PostData(2,fCFManager->GetParticleContainer()) ;
737 //PostData(3,fCorrelation) ;
741 //___________________________________________________________________________
742 void AliCFHeavyFlavourTaskMultiVarMultiStep::Terminate(Option_t*)
744 // The Terminate() function is the last function to be called during
745 // a query. It always runs on the client, it can be used to present
746 // the results graphically or save the results to file.
748 AliAnalysisTaskSE::Terminate();
750 AliInfo(Form("Found %i MC particles that are D0 in MC, in %d events",fCountMC,fEvents));
751 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, in %d events",fCountAcc,fEvents));
752 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fEvents));
753 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));
754 AliInfo(Form("Found %i reco D0 that are decaying in K+pi, in %d events",fCountReco,fEvents));
755 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));
756 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and have at least %d clusters in ITS, in %d events",fCountRecoITSClusters,fMinITSClusters,fEvents));
757 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and satisfy PPR cuts, in %d events",fCountRecoPPR,fEvents));
758 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and satisfy PID cuts, in %d events",fCountRecoPID,fEvents));
760 // draw some example plots....
762 // AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
763 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
765 printf("CONTAINER NOT FOUND\n");
768 // projecting the containers to obtain histograms
769 // first argument = variable, second argument = step
772 TH1D* h00 = cont->ShowProjection(0,0) ; // pt
773 TH1D* h10 = cont->ShowProjection(1,0) ; // rapidity
774 TH1D* h20 = cont->ShowProjection(2,0) ; // cosThetaStar
775 TH1D* h30 = cont->ShowProjection(3,0) ; // pTpi
776 TH1D* h40 = cont->ShowProjection(4,0) ; // pTK
777 TH1D* h50 = cont->ShowProjection(5,0) ; // cT
778 TH1D* h60 = cont->ShowProjection(6,0) ; // dca
779 TH1D* h70 = cont->ShowProjection(7,0) ; // d0pi
780 TH1D* h80 = cont->ShowProjection(8,0) ; // d0K
781 TH1D* h90 = cont->ShowProjection(9,0) ; // d0xd0
782 TH1D* h100 = cont->ShowProjection(10,0) ; // cosPointingAngle
783 TH1D* h110 = cont->ShowProjection(11,0) ; // phi
785 // MC-Acceptance level
786 TH1D* h01 = cont->ShowProjection(0,1) ; // pt
787 TH1D* h11 = cont->ShowProjection(1,1) ; // rapidity
788 TH1D* h21 = cont->ShowProjection(2,1) ; // cosThetaStar
789 TH1D* h31 = cont->ShowProjection(3,1) ; // pTpi
790 TH1D* h41 = cont->ShowProjection(4,1) ; // pTK
791 TH1D* h51 = cont->ShowProjection(5,1) ; // cT
792 TH1D* h61 = cont->ShowProjection(6,1) ; // dca
793 TH1D* h71 = cont->ShowProjection(7,1) ; // d0pi
794 TH1D* h81 = cont->ShowProjection(8,1) ; // d0K
795 TH1D* h91 = cont->ShowProjection(9,1) ; // d0xd0
796 TH1D* h101 = cont->ShowProjection(10,1) ; // cosPointingAngle
797 TH1D* h111 = cont->ShowProjection(11,1) ; // phi
800 TH1D* h04 = cont->ShowProjection(0,4) ; // pt
801 TH1D* h14 = cont->ShowProjection(1,4) ; // rapidity
802 TH1D* h24 = cont->ShowProjection(2,4) ; // cosThetaStar
803 TH1D* h34 = cont->ShowProjection(3,4) ; // pTpi
804 TH1D* h44 = cont->ShowProjection(4,4) ; // pTK
805 TH1D* h54 = cont->ShowProjection(5,4) ; // cT
806 TH1D* h64 = cont->ShowProjection(6,4) ; // dca
807 TH1D* h74 = cont->ShowProjection(7,4) ; // d0pi
808 TH1D* h84 = cont->ShowProjection(8,4) ; // d0K
809 TH1D* h94 = cont->ShowProjection(9,4) ; // d0xd0
810 TH1D* h104 = cont->ShowProjection(10,4) ; // cosPointingAngle
811 TH1D* h114 = cont->ShowProjection(11,4) ; // phi
813 h00->SetTitle("pT_D0 (GeV/c)");
814 h10->SetTitle("rapidity");
815 h20->SetTitle("cosThetaStar");
816 h30->SetTitle("pT_pi (GeV/c)");
817 h40->SetTitle("pT_K (Gev/c)");
818 h50->SetTitle("cT (#mum)");
819 h60->SetTitle("dca (#mum)");
820 h70->SetTitle("d0_pi (#mum)");
821 h80->SetTitle("d0_K (#mum)");
822 h90->SetTitle("d0xd0 (#mum^2)");
823 h100->SetTitle("cosPointingAngle");
824 h100->SetTitle("phi (rad)");
826 h00->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
827 h10->GetXaxis()->SetTitle("rapidity");
828 h20->GetXaxis()->SetTitle("cosThetaStar");
829 h30->GetXaxis()->SetTitle("pT_pi (GeV/c)");
830 h40->GetXaxis()->SetTitle("pT_K (Gev/c)");
831 h50->GetXaxis()->SetTitle("cT (#mum)");
832 h60->GetXaxis()->SetTitle("dca (#mum)");
833 h70->GetXaxis()->SetTitle("d0_pi (#mum)");
834 h80->GetXaxis()->SetTitle("d0_K (#mum)");
835 h90->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
836 h100->GetXaxis()->SetTitle("cosPointingAngle");
837 h110->GetXaxis()->SetTitle("phi (rad)");
839 h01->SetTitle("pT_D0 (GeV/c)");
840 h11->SetTitle("rapidity");
841 h21->SetTitle("cosThetaStar");
842 h31->SetTitle("pT_pi (GeV/c)");
843 h41->SetTitle("pT_K (Gev/c)");
844 h51->SetTitle("cT (#mum)");
845 h61->SetTitle("dca (#mum)");
846 h71->SetTitle("d0_pi (#mum)");
847 h81->SetTitle("d0_K (#mum)");
848 h91->SetTitle("d0xd0 (#mum^2)");
849 h101->SetTitle("cosPointingAngle");
850 h111->GetXaxis()->SetTitle("phi (rad)");
852 h01->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
853 h11->GetXaxis()->SetTitle("rapidity");
854 h21->GetXaxis()->SetTitle("cosThetaStar");
855 h31->GetXaxis()->SetTitle("pT_pi (GeV/c)");
856 h41->GetXaxis()->SetTitle("pT_K (Gev/c)");
857 h51->GetXaxis()->SetTitle("cT (#mum)");
858 h61->GetXaxis()->SetTitle("dca (#mum)");
859 h71->GetXaxis()->SetTitle("d0_pi (#mum)");
860 h81->GetXaxis()->SetTitle("d0_K (#mum)");
861 h91->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
862 h101->GetXaxis()->SetTitle("cosPointingAngle");
863 h111->GetXaxis()->SetTitle("phi (rad)");
865 h04->SetTitle("pT_D0 (GeV/c)");
866 h14->SetTitle("rapidity");
867 h24->SetTitle("cosThetaStar");
868 h34->SetTitle("pT_pi (GeV/c)");
869 h44->SetTitle("pT_K (Gev/c)");
870 h54->SetTitle("cT (#mum)");
871 h64->SetTitle("dca (#mum)");
872 h74->SetTitle("d0_pi (#mum)");
873 h84->SetTitle("d0_K (#mum)");
874 h94->SetTitle("d0xd0 (#mum^2)");
875 h104->SetTitle("cosPointingAngle");
876 h114->GetXaxis()->SetTitle("phi (rad)");
878 h04->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
879 h14->GetXaxis()->SetTitle("rapidity");
880 h24->GetXaxis()->SetTitle("cosThetaStar");
881 h34->GetXaxis()->SetTitle("pT_pi (GeV/c)");
882 h44->GetXaxis()->SetTitle("pT_K (Gev/c)");
883 h54->GetXaxis()->SetTitle("cT (#mum)");
884 h64->GetXaxis()->SetTitle("dca (#mum)");
885 h74->GetXaxis()->SetTitle("d0_pi (#mum)");
886 h84->GetXaxis()->SetTitle("d0_K (#mum)");
887 h94->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
888 h104->GetXaxis()->SetTitle("cosPointingAngle");
889 h114->GetXaxis()->SetTitle("phi (rad)");
891 Double_t max0 = h00->GetMaximum();
892 Double_t max1 = h10->GetMaximum();
893 Double_t max2 = h20->GetMaximum();
894 Double_t max3 = h30->GetMaximum();
895 Double_t max4 = h40->GetMaximum();
896 Double_t max5 = h50->GetMaximum();
897 Double_t max6 = h60->GetMaximum();
898 Double_t max7 = h70->GetMaximum();
899 Double_t max8 = h80->GetMaximum();
900 Double_t max9 = h90->GetMaximum();
901 Double_t max10 = h100->GetMaximum();
902 Double_t max11 = h110->GetMaximum();
904 h00->GetYaxis()->SetRangeUser(0,max0*1.2);
905 h10->GetYaxis()->SetRangeUser(0,max1*1.2);
906 h20->GetYaxis()->SetRangeUser(0,max2*1.2);
907 h30->GetYaxis()->SetRangeUser(0,max3*1.2);
908 h40->GetYaxis()->SetRangeUser(0,max4*1.2);
909 h50->GetYaxis()->SetRangeUser(0,max5*1.2);
910 h60->GetYaxis()->SetRangeUser(0,max6*1.2);
911 h70->GetYaxis()->SetRangeUser(0,max7*1.2);
912 h80->GetYaxis()->SetRangeUser(0,max8*1.2);
913 h90->GetYaxis()->SetRangeUser(0,max9*1.2);
914 h100->GetYaxis()->SetRangeUser(0,max10*1.2);
915 h110->GetYaxis()->SetRangeUser(0,max11*1.2);
917 h01->GetYaxis()->SetRangeUser(0,max0*1.2);
918 h11->GetYaxis()->SetRangeUser(0,max1*1.2);
919 h21->GetYaxis()->SetRangeUser(0,max2*1.2);
920 h31->GetYaxis()->SetRangeUser(0,max3*1.2);
921 h41->GetYaxis()->SetRangeUser(0,max4*1.2);
922 h51->GetYaxis()->SetRangeUser(0,max5*1.2);
923 h61->GetYaxis()->SetRangeUser(0,max6*1.2);
924 h71->GetYaxis()->SetRangeUser(0,max7*1.2);
925 h81->GetYaxis()->SetRangeUser(0,max8*1.2);
926 h91->GetYaxis()->SetRangeUser(0,max9*1.2);
927 h101->GetYaxis()->SetRangeUser(0,max10*1.2);
928 h111->GetYaxis()->SetRangeUser(0,max11*1.2);
930 h00->SetMarkerStyle(20);
931 h10->SetMarkerStyle(24);
932 h20->SetMarkerStyle(21);
933 h30->SetMarkerStyle(25);
934 h40->SetMarkerStyle(27);
935 h50->SetMarkerStyle(28);
936 h60->SetMarkerStyle(20);
937 h70->SetMarkerStyle(24);
938 h80->SetMarkerStyle(21);
939 h90->SetMarkerStyle(25);
940 h100->SetMarkerStyle(27);
941 h110->SetMarkerStyle(28);
943 h00->SetMarkerColor(2);
944 h10->SetMarkerColor(2);
945 h20->SetMarkerColor(2);
946 h30->SetMarkerColor(2);
947 h40->SetMarkerColor(2);
948 h50->SetMarkerColor(2);
949 h60->SetMarkerColor(2);
950 h70->SetMarkerColor(2);
951 h80->SetMarkerColor(2);
952 h90->SetMarkerColor(2);
953 h100->SetMarkerColor(2);
954 h110->SetMarkerColor(2);
956 h01->SetMarkerStyle(20) ;
957 h11->SetMarkerStyle(24) ;
958 h21->SetMarkerStyle(21) ;
959 h31->SetMarkerStyle(25) ;
960 h41->SetMarkerStyle(27) ;
961 h51->SetMarkerStyle(28) ;
962 h61->SetMarkerStyle(20);
963 h71->SetMarkerStyle(24);
964 h81->SetMarkerStyle(21);
965 h91->SetMarkerStyle(25);
966 h101->SetMarkerStyle(27);
967 h111->SetMarkerStyle(28);
969 h01->SetMarkerColor(8);
970 h11->SetMarkerColor(8);
971 h21->SetMarkerColor(8);
972 h31->SetMarkerColor(8);
973 h41->SetMarkerColor(8);
974 h51->SetMarkerColor(8);
975 h61->SetMarkerColor(8);
976 h71->SetMarkerColor(8);
977 h81->SetMarkerColor(8);
978 h91->SetMarkerColor(8);
979 h101->SetMarkerColor(8);
980 h111->SetMarkerColor(8);
982 h04->SetMarkerStyle(20) ;
983 h14->SetMarkerStyle(24) ;
984 h24->SetMarkerStyle(21) ;
985 h34->SetMarkerStyle(25) ;
986 h44->SetMarkerStyle(27) ;
987 h54->SetMarkerStyle(28) ;
988 h64->SetMarkerStyle(20);
989 h74->SetMarkerStyle(24);
990 h84->SetMarkerStyle(21);
991 h94->SetMarkerStyle(25);
992 h104->SetMarkerStyle(27);
993 h114->SetMarkerStyle(28);
995 h04->SetMarkerColor(4);
996 h14->SetMarkerColor(4);
997 h24->SetMarkerColor(4);
998 h34->SetMarkerColor(4);
999 h44->SetMarkerColor(4);
1000 h54->SetMarkerColor(4);
1001 h64->SetMarkerColor(4);
1002 h74->SetMarkerColor(4);
1003 h84->SetMarkerColor(4);
1004 h94->SetMarkerColor(4);
1005 h104->SetMarkerColor(4);
1006 h114->SetMarkerColor(4);
1008 gStyle->SetCanvasColor(0);
1009 gStyle->SetFrameFillColor(0);
1010 gStyle->SetTitleFillColor(0);
1011 gStyle->SetStatColor(0);
1013 // drawing in 2 separate canvas for a matter of clearity
1014 TCanvas * c1 =new TCanvas("c1","pT, rapidiy, cosThetaStar",1100,1600);
1046 TCanvas * c2 =new TCanvas("c2","pTpi, pTK, cT",1100,1600);
1077 TCanvas * c3 =new TCanvas("c3","dca, d0pi, d0K",1100,1600);
1108 TCanvas * c4 =new TCanvas("c4","d0xd0, cosPointingAngle, phi",1100,1600);
1139 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
1141 TH2D* corr1 =hcorr->Projection(0,2);
1142 TH2D* corr2 = hcorr->Projection(1,3);
1144 TCanvas * c7 =new TCanvas("c7","",800,400);
1147 corr1->Draw("text");
1149 corr2->Draw("text");
1152 TString projectionFilename="CFtaskHFprojection";
1153 if(fKeepD0fromBOnly) projectionFilename+="_KeepD0fromBOnly";
1154 else if(fKeepD0fromB) projectionFilename+="_KeepD0fromB";
1155 projectionFilename+=".root";
1156 TFile* fileProjection = new TFile(projectionFilename.Data(),"RECREATE");
1160 h00->Write("pT_D0_step0");
1161 h10->Write("rapidity_step0");
1162 h20->Write("cosThetaStar_step0");
1163 h30->Write("pT_pi_step0");
1164 h40->Write("pT_K_step0");
1165 h50->Write("cT_step0");
1166 h60->Write("dca_step0");
1167 h70->Write("d0_pi_step0");
1168 h80->Write("d0_K_step0");
1169 h90->Write("d0xd0_step0");
1170 h100->Write("cosPointingAngle_step0");
1171 h110->Write("phi_step0");
1173 h01->Write("pT_D0_step1");
1174 h11->Write("rapidity_step1");
1175 h21->Write("cosThetaStar_step1");
1176 h31->Write("pT_pi_step1");
1177 h41->Write("pT_K_step1");
1178 h51->Write("cT_step1");
1179 h61->Write("dca_step1");
1180 h71->Write("d0_pi_step1");
1181 h81->Write("d0_K_step1");
1182 h91->Write("d0xd0_step1");
1183 h101->Write("cosPointingAngle_step1");
1184 h111->Write("phi_step1");
1186 h04->Write("pT_D0_step2");
1187 h14->Write("rapidity_step2");
1188 h24->Write("cosThetaStar_step2");
1189 h34->Write("pT_pi_step2");
1190 h44->Write("pT_K_step2");
1191 h54->Write("cT_step2");
1192 h64->Write("dca_step2");
1193 h74->Write("d0_pi_step2");
1194 h80->Write("d0_K_step2");
1195 h94->Write("d0xd0_step2");
1196 h104->Write("cosPointingAngle_step2");
1197 h114->Write("phi_step2");
1199 fileProjection->Close();
1204 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
1205 c2->SaveAs("Plots/pTpi_pTK_cT.eps");
1206 c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
1207 c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
1209 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
1210 c2->SaveAs("Plots/pTpi_pTK_cT.gif");
1211 c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
1212 c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
1216 //___________________________________________________________________________
1217 void AliCFHeavyFlavourTaskMultiVarMultiStep::UserCreateOutputObjects() {
1218 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1219 //TO BE SET BEFORE THE EXECUTION OF THE TASK
1221 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1225 fHistEventsProcessed = new TH1I("CFHFchist0","",1,0,1) ;
1228 //___________________________________________________________________________
1229 Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CosThetaStar(const AliAODMCParticle* mcPart, const AliAODMCParticle* mcPartDaughter0, const AliAODMCParticle* mcPartDaughter1) const {
1232 // to calculate cos(ThetaStar) for generated particle
1233 // using the K, since mcPartDaughter0 and mcPartDaughter1 always correspond to K and pi respectively
1234 // (see where the function is called)
1237 Int_t pdgvtx = mcPart->GetPdgCode();
1238 /* if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1239 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1240 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1241 AliInfo(Form("D0, with pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1242 AliDebug(2,"This is a D0");
1243 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1244 mcPartDaughter0 = mcPartDaughter1;
1245 mcPartDaughter1 = mcPartdummy;
1251 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1252 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1253 if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1257 AliDebug(2,"D0bar");
1259 if (pdgprong0 == 211){
1260 AliDebug(2,Form("pdgprong0 = %d, pdgprong1 = %d, switching...",pdgprong0,pdgprong1));
1261 const AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1262 mcPartDaughter0 = mcPartDaughter1;
1263 mcPartDaughter1 = mcPartdummy;
1264 pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1265 pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1268 AliDebug(2,Form("After checking, pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1269 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdgvtx))->Mass();
1271 massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
1272 massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
1274 Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)-4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx);
1276 Double_t px = mcPartDaughter0->Px()+mcPartDaughter1->Px();
1277 Double_t py = mcPartDaughter0->Py()+mcPartDaughter1->Py();
1278 Double_t pz = mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
1279 Double_t p = TMath::Sqrt(px*px+py*py+pz*pz);
1280 Double_t e = TMath::Sqrt(massvtx*massvtx+p*p);
1281 Double_t beta = p/e;
1282 Double_t gamma = e/massvtx;
1283 TVector3 mom(mcPartDaughter0->Px(),mcPartDaughter0->Py(),mcPartDaughter0->Pz());
1284 TVector3 momTot(mcPartDaughter0->Px()+mcPartDaughter1->Px(),mcPartDaughter0->Py()+mcPartDaughter1->Py(),mcPartDaughter0->Pz()+mcPartDaughter1->Pz());
1286 Double_t qlprong = mom.Dot(momTot)/momTot.Mag(); // analog to AliAODRecoDecay::QlProng(0)
1288 AliDebug(2,Form("pStar = %f, beta = %f, gamma = %f, qlprong = %f, massp[0] = %f", pStar, beta, gamma, qlprong, massp[0]));
1289 Double_t cts = (qlprong/gamma-beta*TMath::Sqrt(pStar*pStar+massp[0]*massp[0]))/pStar;
1290 AliDebug(2,Form("cts = %f", cts));
1293 //___________________________________________________________________________
1294 Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CT(const AliAODMCParticle* mcPart, const AliAODMCParticle* mcPartDaughter0, const AliAODMCParticle* mcPartDaughter1) const {
1297 // to calculate cT for generated particle
1300 Double_t xmcPart[3] = {0,0,0};
1301 Double_t xdaughter0[3] = {0,0,0};
1302 Double_t xdaughter1[3] = {0,0,0};
1303 mcPart->XvYvZv(xmcPart); // cm
1304 mcPartDaughter0->XvYvZv(xdaughter0); // cm
1305 mcPartDaughter1->XvYvZv(xdaughter1); //cm
1306 Double_t prodVtxD0 = TMath::Sqrt(xmcPart[0]*xmcPart[0]+
1307 xmcPart[1]*xmcPart[1]+
1308 xmcPart[2]*xmcPart[2]);
1309 Double_t prodVtxDaughter0 = TMath::Sqrt(xdaughter0[0]*xdaughter0[0]+
1310 xdaughter0[1]*xdaughter0[1]+
1311 xdaughter0[2]*xdaughter0[2]);
1312 Double_t prodVtxDaughter1 = TMath::Sqrt(xdaughter1[0]*xdaughter1[0]+
1313 xdaughter1[1]*xdaughter1[1]+
1314 xdaughter1[2]*xdaughter1[2]);
1316 AliDebug(2, Form("D0: x = %f, y = %f, z = %f, production vertex distance = %f (cm), %f (micron)", xmcPart[0], xmcPart[1], xmcPart[2], prodVtxD0, prodVtxD0*1.E4));
1317 AliDebug(2, Form("Daughter0: x = %f, y = %f, z = %f, production vertex distance = %f (cm) %f (micron)", xdaughter0[0], xdaughter0[1], xdaughter0[2], prodVtxDaughter0, prodVtxDaughter0*1E4));
1318 AliDebug(2, Form("Daughter1: x = %f, y = %f, z = %f, production vertex distance = %f (cm) %f (micron)", xdaughter1[0], xdaughter1[1], xdaughter1[2], prodVtxDaughter1, prodVtxDaughter1*1.E4));
1320 Double_t cT0 = TMath::Sqrt((xdaughter0[0]-xmcPart[0])*(xdaughter0[0]-xmcPart[0])+
1321 (xdaughter0[1]-xmcPart[1])*(xdaughter0[1]-xmcPart[1])+
1322 (xdaughter0[2]-xmcPart[2])*(xdaughter0[2]-xmcPart[2]));
1324 Double_t cT1 = TMath::Sqrt((xdaughter1[0]-xmcPart[0])*(xdaughter1[0]-xmcPart[0])+
1325 (xdaughter1[1]-xmcPart[1])*(xdaughter1[1]-xmcPart[1])+
1326 (xdaughter1[2]-xmcPart[2])*(xdaughter1[2]-xmcPart[2]));
1328 if ((cT0 - cT1)>1E-5) {
1329 AliWarning(Form("cT from daughter 0 (%f) different from cT from daughter 1 (%f)! Using cT from daughter 0, but PLEASE, CHECK....",cT0,cT1));
1332 // calculating cT from cT0
1334 Double_t mass = TDatabasePDG::Instance()->GetParticle(mcPart->GetPdgCode())->Mass(); // mcPart->GetPdgCode() should return +/- 421 for the D0/D0bar
1335 Double_t p = mcPart-> P();
1336 Double_t cT = cT0*mass/p;
1337 AliDebug(2, Form("cT from daughter 0 = %f (micron)", cT0*1E4));
1338 AliDebug(2, Form("cT from daughter 1 = %f (micron)", cT1*1E4));
1339 AliDebug(2, Form("cT (with mass = %f and p = %f) = %f (micron)", mass, p, cT*1E4));
1342 //________________________________________________________________________________
1343 Bool_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, const TClonesArray* mcArray, Double_t* vectorMC) const {
1346 // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
1349 Bool_t isOk = kFALSE;
1351 // check whether the D0 decays in pi+K
1352 // to be added!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1353 // could use a cut in the CF, but not clear how to define a TDecayChannel
1354 // implemented for the time being as a cut on the number of daughters - see later when
1355 // getting the daughters
1357 // getting the daughters
1358 Int_t daughter0 = mcPart->GetDaughter(0);
1359 Int_t daughter1 = mcPart->GetDaughter(1);
1360 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
1361 if (daughter0 == 0 || daughter1 == 0) {
1362 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
1365 if (TMath::Abs(daughter1 - daughter0) != 1) {
1366 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
1369 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
1370 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
1371 if (!mcPartDaughter0 || !mcPartDaughter1) {
1372 AliWarning("At least one Daughter Particle not found in tree, skipping");
1375 if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==321 &&
1376 TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
1377 !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
1378 TMath::Abs(mcPartDaughter1->GetPdgCode())==321)) {
1379 AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
1383 Double_t vtx1[3] = {0,0,0}; // primary vertex
1384 Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
1385 Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
1386 mcPart->XvYvZv(vtx1); // cm
1387 // getting vertex from daughters
1388 mcPartDaughter0->XvYvZv(vtx2daughter0); // cm
1389 mcPartDaughter1->XvYvZv(vtx2daughter1); //cm
1390 if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
1391 AliError("Daughters have different secondary vertex, skipping the track");
1396 // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
1397 AliAODMCParticle* positiveDaugh = mcPartDaughter0;
1398 AliAODMCParticle* negativeDaugh = mcPartDaughter1;
1399 if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
1400 // inverting in case the positive daughter is the second one
1401 positiveDaugh = mcPartDaughter1;
1402 negativeDaugh = mcPartDaughter0;
1404 // getting the momentum from the daughters
1405 Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
1406 Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
1407 Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
1409 Double_t d0[2] = {0.,0.};
1411 AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
1413 Double_t cosThetaStar = 0.;
1414 Double_t cosThetaStarD0 = 0.;
1415 Double_t cosThetaStarD0bar = 0.;
1416 cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
1417 cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
1418 if (mcPart->GetPdgCode() == 421){ // D0
1419 AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1420 cosThetaStar = cosThetaStarD0;
1422 else if (mcPart->GetPdgCode() == -421){ // D0bar{
1423 AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1424 cosThetaStar = cosThetaStarD0bar;
1427 AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
1430 if (cosThetaStar < -1 || cosThetaStar > 1) {
1431 AliWarning("Invalid value for cosine Theta star, returning");
1435 // calculate cos(Theta*) according to the method implemented herein
1437 Double_t cts = 9999.;
1438 cts = CosThetaStar(mcPart, mcPartDaughter0, mcPartDaughter1);
1439 if (cts < -1 || cts > 1) {
1440 AliWarning("Invalid value for cosine Theta star from AliCFHeavyFlavourTaskMultiVarMultiStep method");
1442 if (TMath::Abs(cts - cosThetaStar)>0.001) {
1443 AliError(Form("cosThetaStar automatically calculated different from that manually calculated!!! cosThetaStar = %f, cosThetaStar = %f", cosThetaStar,cts));
1446 Double_t cT = decay->Ct(421);
1448 // calculate cT from the method implemented herein
1450 cT1 = CT(mcPart, mcPartDaughter0, mcPartDaughter1);
1452 if (TMath::Abs(cT1 - cT)>0.001) {
1453 AliError(Form("cT automatically calculated different from that manually calculated!!! cT = %f, cT1 = %f",cT,cT1));
1456 // get the pT of the daughters
1461 if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
1462 pTpi = mcPartDaughter0->Pt();
1463 pTK = mcPartDaughter1->Pt();
1466 pTpi = mcPartDaughter1->Pt();
1467 pTK = mcPartDaughter0->Pt();
1470 vectorMC[0] = mcPart->Pt();
1471 vectorMC[1] = mcPart->Y() ;
1472 vectorMC[2] = cosThetaStar ;
1473 vectorMC[3] = pTpi ;
1475 vectorMC[5] = cT*1.E4 ; // in micron
1476 vectorMC[6] = mcPart->Phi() ;
1480 //_________________________________________________________________________________________________
1481 Int_t AliCFHeavyFlavourTaskMultiVarMultiStep::CheckOrigin(const AliAODMCParticle* mcPart, const TClonesArray* mcArray)const{
1484 // checking whether the very mother of the D0 is a charm or a bottom quark
1487 Int_t pdgGranma = 0;
1489 mother = mcPart->GetMother();
1493 AliDebug(2,Form("mother at step %d = %d", istep, mother));
1494 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
1495 if(!mcGranma) break;
1496 pdgGranma = mcGranma->GetPdgCode();
1497 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1498 Int_t abspdgGranma = TMath::Abs(pdgGranma);
1499 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1502 mother = mcGranma->GetMother();
1506 //__________________________________________________________________________________________________
1507 Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetWeight(Float_t pt){
1510 // calculating the weight to fill the container
1514 // p0 = 1.63297e-01 --> 0.322643
1520 // p0 = 1.85906e-01 --> 0.36609
1525 Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
1526 Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
1528 Double_t dndptFunc1 = DodNdptFit(pt,func1);
1529 Double_t dndptFunc2 = DodNdptFit(pt,func2);
1530 AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndptFunc1,dndptFunc2,dndptFunc1/dndptFunc2));
1531 return dndptFunc1/dndptFunc2;
1534 //__________________________________________________________________________________________________
1535 Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::DodNdptFit(Float_t pt, const Double_t* par)const{
1538 // calculating dNdpt
1541 Double_t denom = TMath::Power((pt/par[1]), par[3] );
1542 Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);