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 // 13 variables used: pt, y, cosThetaStar, ptPi, ptK, ct,
21 // dca, d0Pi, d0K, d0Pixd0K, cosPointingAngle, phi, z
23 //-----------------------------------------------------------------------
24 // Author : C. Zampolli, CERN
25 //-----------------------------------------------------------------------
26 //-----------------------------------------------------------------------
27 // Base class for HF Unfolding (pt and eta)
28 // correlation matrix filled at Acceptance and PPR level
29 // Author: A.Grelli , Utrecht - agrelli@uu.nl
30 //-----------------------------------------------------------------------
32 #include <TParticle.h>
33 #include <TDatabasePDG.h>
38 #include "AliCFHeavyFlavourTaskMultiVarMultiStep.h"
40 #include "AliMCEvent.h"
41 #include "AliCFManager.h"
42 #include "AliCFContainer.h"
44 #include "AliAnalysisManager.h"
45 #include "AliAODHandler.h"
46 #include "AliAODEvent.h"
47 #include "AliAODRecoDecay.h"
48 #include "AliAODRecoDecayHF.h"
49 #include "AliAODRecoDecayHF2Prong.h"
50 #include "AliAODMCParticle.h"
51 #include "AliAODMCHeader.h"
52 #include "AliESDtrack.h"
53 #include "AliRDHFCutsD0toKpi.h"
55 #include "THnSparse.h"
57 #include "AliAnalysisDataSlot.h"
58 #include "AliAnalysisDataContainer.h"
60 //__________________________________________________________________________
61 AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep() :
65 fHistEventsProcessed(0x0),
73 fCountRecoITSClusters(0),
77 fFillFromGenerated(kFALSE),
79 fAcceptanceUnf(kTRUE),
81 fKeepD0fromBOnly(kFALSE),
90 //___________________________________________________________________________
91 AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(const Char_t* name, AliRDHFCutsD0toKpi* cuts) :
92 AliAnalysisTaskSE(name),
95 fHistEventsProcessed(0x0),
103 fCountRecoITSClusters(0),
107 fFillFromGenerated(kFALSE),
109 fAcceptanceUnf(kTRUE),
110 fKeepD0fromB(kFALSE),
111 fKeepD0fromBOnly(kFALSE),
117 // Constructor. Initialization of Inputs and Outputs
119 Info("AliCFHeavyFlavourTaskMultiVarMultiStep","Calling Constructor");
121 DefineInput(0) and DefineOutput(0)
122 are taken care of by AliAnalysisTaskSE constructor
124 DefineOutput(1,TH1I::Class());
125 DefineOutput(2,AliCFContainer::Class());
126 DefineOutput(3,THnSparseD::Class());
127 DefineOutput(4,AliRDHFCutsD0toKpi::Class());
132 //___________________________________________________________________________
133 AliCFHeavyFlavourTaskMultiVarMultiStep& AliCFHeavyFlavourTaskMultiVarMultiStep::operator=(const AliCFHeavyFlavourTaskMultiVarMultiStep& c)
136 // Assignment operator
139 AliAnalysisTaskSE::operator=(c) ;
141 fCFManager = c.fCFManager;
142 fHistEventsProcessed = c.fHistEventsProcessed;
148 //___________________________________________________________________________
149 AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(const AliCFHeavyFlavourTaskMultiVarMultiStep& c) :
150 AliAnalysisTaskSE(c),
152 fCFManager(c.fCFManager),
153 fHistEventsProcessed(c.fHistEventsProcessed),
154 fCorrelation(c.fCorrelation),
155 fCountMC(c.fCountMC),
156 fCountAcc(c.fCountAcc),
157 fCountVertex(c.fCountVertex),
158 fCountRefit(c.fCountRefit),
159 fCountReco(c.fCountReco),
160 fCountRecoAcc(c.fCountRecoAcc),
161 fCountRecoITSClusters(c.fCountRecoITSClusters),
162 fCountRecoPPR(c.fCountRecoPPR),
163 fCountRecoPID(c.fCountRecoPID),
165 fFillFromGenerated(c.fFillFromGenerated),
166 fMinITSClusters(c.fMinITSClusters),
167 fAcceptanceUnf(c.fAcceptanceUnf),
168 fKeepD0fromB(c.fKeepD0fromB),
169 fKeepD0fromBOnly(c.fKeepD0fromBOnly),
171 fUseWeight(c.fUseWeight),
179 //___________________________________________________________________________
180 AliCFHeavyFlavourTaskMultiVarMultiStep::~AliCFHeavyFlavourTaskMultiVarMultiStep() {
184 if (fCFManager) delete fCFManager ;
185 if (fHistEventsProcessed) delete fHistEventsProcessed ;
186 if (fCorrelation) delete fCorrelation ;
187 // if (fCuts) delete fCuts ;
190 //________________________________________________________________________
191 void AliCFHeavyFlavourTaskMultiVarMultiStep::Init(){
197 if(fDebug > 1) printf("AliCFHeavyFlavourTaskMultiVarMultiStep::Init() \n");
199 AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*fCuts);
200 const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
201 copyfCuts->SetName(nameoutput);
203 PostData(4,copyfCuts);
208 //_________________________________________________
209 void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
212 // Main loop function
215 PostData(1,fHistEventsProcessed) ;
216 PostData(2,fCFManager->GetParticleContainer()) ;
217 PostData(3,fCorrelation) ;
219 AliESDtrackCuts* trackCuts = fCuts->GetTrackCuts(); // track cuts
221 if (fFillFromGenerated){
222 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
226 Error("UserExec","NO EVENT FOUND!");
230 // check that the fKeepD0fromB flag is set to true when the fKeepD0fromBOnly flag is true
231 if(fKeepD0fromBOnly) {
233 if(fEvents<2) AliInfo(Form("Both fKeepD0fromB and fKeepD0fromBOnly flags are true, looking _ONLY_ at D0 FROM B"));
236 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
238 TClonesArray *arrayD0toKpi=0;
240 if(!aodEvent && AODEvent() && IsStandardAOD()) {
241 // In case there is an AOD handler writing a standard AOD, use the AOD
242 // event in memory rather than the input (ESD) event.
243 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
244 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
245 // have to taken from the AOD event hold by the AliAODExtension
246 AliAODHandler* aodHandler = (AliAODHandler*)
247 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
248 if(aodHandler->GetExtensions()) {
249 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
250 AliAODEvent *aodFromExt = ext->GetAOD();
251 arrayD0toKpi=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
254 arrayD0toKpi=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
259 AliError("Could not find array of HF vertices");
263 // fix for temporary bug in ESDfilter
264 // the AODs with null vertex pointer didn't pass the PhysSel
265 if(!aodEvent->GetPrimaryVertex()) return;
268 if (fEvents%10000 ==0) AliDebug(2,Form("Event %d",fEvents));
270 fCFManager->SetRecEventInfo(aodEvent);
271 fCFManager->SetMCEventInfo(aodEvent);
273 // MC-event selection
274 Double_t containerInput[13] ;
275 Double_t containerInputMC[13] ;
277 //loop on the MC event
279 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
281 AliError("Could not find Monte-Carlo in AOD");
286 Int_t icountReco = 0;
287 Int_t icountVertex = 0;
288 Int_t icountRefit = 0;
289 Int_t icountRecoAcc = 0;
290 Int_t icountRecoITSClusters = 0;
291 Int_t icountRecoPPR = 0;
292 Int_t icountRecoPID = 0;
294 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
296 AliError("Could not find MC Header in AOD");
302 // AOD primary vertex
303 AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
305 AliError("There is no primary vertex !");
308 Double_t zPrimVertex = vtx1->GetZ();
309 Double_t zMCVertex = mcHeader->GetVtxZ();
310 Bool_t vtxFlag = kTRUE;
311 TString title=vtx1->GetTitle();
312 if(!title.Contains("VertexerTracks")) vtxFlag=kFALSE;
314 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
315 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
316 if (mcPart->GetPdgCode() == 4) cquarks++;
317 if (mcPart->GetPdgCode() == -4) cquarks++;
319 AliWarning("Particle not found in tree, skipping");
323 // check the MC-level cuts
325 if (!fCFManager->CheckParticleCuts(0, mcPart)) continue; // 0 stands for MC level
326 Int_t pdgGranma = CheckOrigin(mcPart, mcArray);
327 Int_t abspdgGranma = TMath::Abs(pdgGranma);
328 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
329 AliDebug(2,Form("Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
330 if (!fKeepD0fromB) continue; // skipping particles that don't come from c quark
332 else { if(fKeepD0fromBOnly) continue; } // skipping particles that don't come from b quark
334 // if (TMath::Abs(pdgGranma)!=4) {
336 // fill the container for Gen-level selection
337 Double_t vectorMC[7] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.};
338 if (GetGeneratedValuesFromMCParticle(mcPart, mcArray, vectorMC)){
339 containerInputMC[0] = vectorMC[0];
340 containerInputMC[1] = vectorMC[1] ;
341 containerInputMC[2] = vectorMC[2] ;
342 containerInputMC[3] = vectorMC[3] ;
343 containerInputMC[4] = vectorMC[4] ;
344 containerInputMC[5] = vectorMC[5] ; // in micron
345 containerInputMC[6] = 0.; // dummy value, meaningless in MC, in micron
346 containerInputMC[7] = 0.; // dummy value, meaningless in MC, in micron
347 containerInputMC[8] = 0.; // dummy value, meaningless in MC, in micron
348 containerInputMC[9] = -100000.; // dummy value, meaningless in MC, in micron^2
349 containerInputMC[10] = 1.01; // dummy value, meaningless in MC
350 containerInputMC[11] = vectorMC[6]; // dummy value, meaningless in MC
351 containerInputMC[12] = zMCVertex; // z of reconstructed of primary vertex
352 if (fUseWeight) fWeight = GetWeight(vectorMC[0]); // setting the weight according to the function defined in AliCFHeavyFlavourTaskMultiVarMultiStep::GetWeight(Float_t pt)
353 AliDebug(3,Form("weight = %f",fWeight));
354 if (!fCuts->IsInFiducialAcceptance(vectorMC[0],vectorMC[1])) continue;
355 if (TMath::Abs(vectorMC[1]) < 0.5) {
356 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc,fWeight);
358 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGenerated,fWeight);
361 // check the MC-Acceptance level cuts
362 // since standard CF functions are not applicable, using Kine Cuts on daughters
364 Int_t daughter0 = mcPart->GetDaughter(0);
365 Int_t daughter1 = mcPart->GetDaughter(1);
366 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
367 if (daughter0 == 0 || daughter1 == 0) {
368 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!! But it should have, this check was already done...");
370 if (TMath::Abs(daughter1 - daughter0) != 1) {
371 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, but it should be, this check was already done...");
373 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
374 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
375 if (!mcPartDaughter0 || !mcPartDaughter1) {
376 AliWarning("At least one Daughter Particle not found in tree, but it should be, this check was already done...");
379 Double_t eta0 = mcPartDaughter0->Eta();
380 Double_t eta1 = mcPartDaughter1->Eta();
381 Double_t y0 = mcPartDaughter0->Y();
382 Double_t y1 = mcPartDaughter1->Y();
383 Double_t pt0 = mcPartDaughter0->Pt();
384 Double_t pt1 = mcPartDaughter1->Pt();
385 AliDebug(2, Form("Daughter 0: eta = %f, y = %f, pt = %f", eta0, y0, pt0));
386 AliDebug(2, Form("Daughter 1: eta = %f, y = %f, pt = %f", eta1, y1, pt1));
387 Bool_t daught0inAcceptance = (TMath::Abs(eta0) < 0.9 && pt0 > 0.1);
388 Bool_t daught1inAcceptance = (TMath::Abs(eta1) < 0.9 && pt1 > 0.1);
389 if (daught0inAcceptance && daught1inAcceptance) {
390 // checking whether the cuts implemented in the CF are equivalent - simply a cross-check
391 AliDebug(2, "Daughter particles in acceptance");
392 if (!fCFManager->CheckParticleCuts(1, mcPartDaughter0)) {
393 AliDebug(2,"Inconsistency with CF cut for daughter 0!");
395 if (!fCFManager->CheckParticleCuts(1, mcPartDaughter1)) {
396 AliDebug(2,"Inconsistency with CF cut for daughter 1!");
398 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance,fWeight);
401 // check on the vertex
402 if (fCuts->IsEventSelected(aodEvent)){
403 AliDebug(2,"Vertex cut passed\n");
404 // filling the container if the vertex is ok
405 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex,fWeight) ;
407 // check on the kTPCrefit and kITSrefit conditions of the daughters
408 Bool_t refitFlag = kTRUE;
409 if (trackCuts->GetRequireTPCRefit() || trackCuts->GetRequireITSRefit()){
410 Int_t foundDaughters = 0;
411 for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
412 AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
413 if ((track->GetLabel() == daughter0) || (track->GetLabel() == daughter1)) {
415 if (trackCuts->GetRequireTPCRefit()) {
416 if(!(track->GetStatus()&AliESDtrack::kTPCrefit)){
421 if (trackCuts->GetRequireITSRefit()) {
422 if(!(track->GetStatus()&AliESDtrack::kITSrefit)){
428 if (foundDaughters == 2){ // both daughters have been checked
432 if (foundDaughters != 2) refitFlag = kFALSE;
435 AliDebug(3,"Refit cut passed\n");
436 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit,fWeight);
440 AliDebug(3,"Refit cut not passed\n");
444 AliDebug(3,"Vertex cut not passed\n");
448 AliDebug(3,"Acceptance cut not passed\n");
452 AliDebug(3,"Problems in filling the container");
457 if (cquarks<2) AliDebug(2,Form("Event found with %d c-quarks", cquarks));
459 AliDebug(2,Form("Found %i MC particles that are D0!!",icountMC));
460 AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Acc cuts!!",icountAcc));
461 AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Vertex cuts!!",icountVertex));
462 AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Refit cuts!!",icountRefit));
464 // Now go to rec level
465 fCountMC += icountMC;
466 fCountAcc += icountAcc;
468 AliDebug(2, Form("Found %d vertices",arrayD0toKpi->GetEntriesFast()));
470 Int_t pdgDgD0toKpi[2]={321,211};
471 Int_t isD0D0bar=1;// 1 for D0, 2 for D0bar, to be used for the PPR and PID selection steps
473 for (Int_t iD0toKpi = 0; iD0toKpi<arrayD0toKpi->GetEntriesFast(); iD0toKpi++) {
475 AliAODRecoDecayHF2Prong* d0tokpi = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->At(iD0toKpi);
476 if(!d0tokpi) continue;
477 Bool_t unsetvtx=kFALSE;
478 if(!d0tokpi->GetOwnPrimaryVtx()) {
479 d0tokpi->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
483 // find associated MC particle
484 Int_t mcLabel = d0tokpi->MatchToMC(421,mcArray,2,pdgDgD0toKpi) ;
487 AliDebug(2,"No MC particle found");
491 AliAODMCParticle* mcVtxHF = (AliAODMCParticle*)mcArray->At(mcLabel);
493 AliWarning("Could not find associated MC in AOD MC tree");
496 if(mcVtxHF->GetPdgCode()==421)isD0D0bar=1;
497 else if(mcVtxHF->GetPdgCode()==-421)isD0D0bar=2;
499 // check whether the daughters have kTPCrefit and kITSrefit set
500 AliAODTrack *track0 = (AliAODTrack*)d0tokpi->GetDaughter(0);
501 AliAODTrack *track1 = (AliAODTrack*)d0tokpi->GetDaughter(1);
502 if( !track0 || !track1 ) {
503 AliWarning("Could not find associated MC daughter tracks in AOD MC tree");
506 if ((trackCuts->GetRequireTPCRefit() && (!(track0->GetStatus()&AliESDtrack::kTPCrefit) || !(track1->GetStatus()&AliESDtrack::kTPCrefit))) ||
507 (trackCuts->GetRequireITSRefit() && (!(track0->GetStatus()&AliESDtrack::kITSrefit) || !(track1->GetStatus()&AliESDtrack::kITSrefit)))){
508 // skipping if at least one daughter does not have kTPCrefit or kITSrefit, if they were required
512 // check on the vertex -- could be moved outside the loop on the reconstructed D0...
513 if(!fCuts->IsEventSelected(aodEvent)) {
514 // skipping cause vertex is not ok
518 const Double_t d0tokpiCuts[9] = {0.3,999999.,1.1,0.,0.,999999.,999999.,999999.,0.};
520 if (!(d0tokpi->SelectD0(&d0tokpiCuts[0],okD0,okD0bar))){
521 // skipping candidate
525 // check if associated MC v0 passes the cuts
526 if (!fCFManager->CheckParticleCuts(0 ,mcVtxHF)) { // 0 stands for MC
527 AliDebug(2, "Skipping the particles due to cuts");
530 Int_t pdgGranma = CheckOrigin(mcVtxHF, mcArray);
531 Int_t abspdgGranma = TMath::Abs(pdgGranma);
532 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
533 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));
534 if (!fKeepD0fromB) continue; // skipping particles that don't come from c quark
536 else { if(fKeepD0fromBOnly) continue; } // skipping particles that don't come from b quark
538 // fill the container...
539 Double_t pt = d0tokpi->Pt();
540 Double_t rapidity = d0tokpi->YD0();
542 Double_t cosThetaStar = 9999.;
545 Double_t dca = d0tokpi->GetDCA();
548 Double_t d0xd0 = d0tokpi->Prodd0d0();
549 Double_t cosPointingAngle = d0tokpi->CosPointingAngle();
550 Double_t phi = d0tokpi->Phi();
551 Int_t pdgCode = mcVtxHF->GetPdgCode();
553 cosThetaStar = d0tokpi->CosThetaStarD0();
554 pTpi = d0tokpi->PtProng(0);
555 pTK = d0tokpi->PtProng(1);
556 d0pi = d0tokpi->Getd0Prong(0);
557 d0K = d0tokpi->Getd0Prong(1);
558 invMass=d0tokpi->InvMassD0();
561 cosThetaStar = d0tokpi->CosThetaStarD0bar();
562 pTpi = d0tokpi->PtProng(1);
563 pTK = d0tokpi->PtProng(0);
564 d0pi = d0tokpi->Getd0Prong(1);
565 d0K = d0tokpi->Getd0Prong(0);
566 invMass=d0tokpi->InvMassD0bar();
569 Double_t cT = d0tokpi->CtD0();
571 if (!fFillFromGenerated){
572 // ...either with reconstructed values....
573 containerInput[0] = pt;
574 containerInput[1] = rapidity;
575 containerInput[2] = cosThetaStar;
576 containerInput[3] = pTpi;
577 containerInput[4] = pTK;
578 containerInput[5] = cT*1.E4; // in micron
579 containerInput[6] = dca*1.E4; // in micron
580 containerInput[7] = d0pi*1.E4; // in micron
581 containerInput[8] = d0K*1.E4; // in micron
582 containerInput[9] = d0xd0*1.E8; // in micron^2
583 containerInput[10] = cosPointingAngle; // in micron
584 containerInput[11] = phi;
585 containerInput[12] = zPrimVertex; // z of reconstructed of primary vertex
588 // ... or with generated values
589 Double_t vectorMC[7] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.};
590 if (GetGeneratedValuesFromMCParticle(mcVtxHF, mcArray, vectorMC)){
591 containerInput[0] = vectorMC[0];
592 containerInput[1] = vectorMC[1] ;
593 containerInput[2] = vectorMC[2] ;
594 containerInput[3] = vectorMC[3] ;
595 containerInput[4] = vectorMC[4] ;
596 containerInput[5] = vectorMC[5] ; // in micron
597 containerInput[6] = 0.; // dummy value, meaningless in MC, in micron
598 containerInput[7] = 0.; // dummy value, meaningless in MC, in micron
599 containerInput[8] = 0.; // dummy value, meaningless in MC, in micron
600 containerInput[9] = 100000.; // dummy value, meaningless in MC, in micron^2
601 containerInput[10] = 1.01; // dummy value, meaningless in MC
602 containerInput[11] = vectorMC[6];
603 containerInput[12] = zMCVertex; // z of reconstructed of primary vertex
606 AliDebug(3,"Problems in filling the container");
611 if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue; // fiducial region
613 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]));
615 AliDebug(2,Form("%d: filling RECO step\n",iD0toKpi));
616 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed,fWeight) ;
619 Float_t etaCutMin, etaCutMax, ptCutMin, ptCutMax;
620 trackCuts->GetEtaRange(etaCutMin,etaCutMax);
621 trackCuts->GetPtRange(ptCutMin,ptCutMax);
622 Bool_t acceptanceProng0 = (d0tokpi->EtaProng(0)>etaCutMin && d0tokpi->EtaProng(0)<etaCutMax && d0tokpi->PtProng(0) > ptCutMin && d0tokpi->PtProng(0) < ptCutMax);
623 Bool_t acceptanceProng1 = (d0tokpi->EtaProng(1)>etaCutMin && d0tokpi->EtaProng(1)<etaCutMax && d0tokpi->PtProng(1) > ptCutMin && d0tokpi->PtProng(1) < ptCutMax);
624 if (acceptanceProng0 && acceptanceProng1) {
625 AliDebug(2,"D0 reco daughters in acceptance");
626 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance,fWeight) ;
631 Double_t fill[4]; //fill response matrix
633 // dimensions 0&1 : pt,eta (Rec)
638 // dimensions 2&3 : pt,eta (MC)
640 fill[2] = mcVtxHF->Pt();
641 fill[3] = mcVtxHF->Y();
643 fCorrelation->Fill(fill);
647 // cut on the min n. of clusters in ITS
648 if (fCuts->IsSelected(d0tokpi,AliRDHFCuts::kTracks)){
649 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
650 icountRecoITSClusters++;
651 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));
653 // setting the use of the PID cut when applying the selection to FALSE - whatever it was. Keeping track of the original value
654 Bool_t iscutusingpid=fCuts->GetIsUsePID();
655 Int_t isselcuts=-1,isselpid=-1;
656 fCuts->SetUsePID(kFALSE);
657 isselcuts = fCuts->IsSelected(d0tokpi,AliRDHFCuts::kCandidate,aodEvent);
658 fCuts->SetUsePID(iscutusingpid); // restoring usage of the PID from the cuts object
659 if (isselcuts == 3 || isselcuts == isD0D0bar){
660 AliDebug(2,"Particle passed PPR cuts (actually cuts for D0 analysis!)");
661 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPPR,fWeight) ;
664 if(!fAcceptanceUnf){ // unfolding
666 Double_t fill[4]; //fill response matrix
668 // dimensions 0&1 : pt,eta (Rec)
673 // dimensions 2&3 : pt,eta (MC)
675 fill[2] = mcVtxHF->Pt();
676 fill[3] = mcVtxHF->Y();
678 fCorrelation->Fill(fill);
682 isselpid = fCuts->IsSelected(d0tokpi,AliRDHFCuts::kPID);
683 if((fCuts->CombineSelectionLevels(3,isselcuts,isselpid)==isD0D0bar)||(fCuts->CombineSelectionLevels(3,isselcuts,isselpid)==3)){
684 AliDebug(2,"Particle passed PID cuts");
685 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPID,fWeight) ;
692 if(unsetvtx) d0tokpi->UnsetOwnPrimaryVtx();
693 } // end loop on D0->Kpi
695 AliDebug(2, Form("Found %i Reco particles that are D0!!",icountReco));
697 fCountReco+= icountReco;
698 fCountVertex+= icountVertex;
699 fCountRefit+= icountRefit;
700 fCountRecoAcc+= icountRecoAcc;
701 fCountRecoITSClusters+= icountRecoITSClusters;
702 fCountRecoPPR+= icountRecoPPR;
703 fCountRecoPID+= icountRecoPID;
705 fHistEventsProcessed->Fill(0);
706 /* PostData(0) is taken care of by AliAnalysisTaskSE */
707 //PostData(1,fHistEventsProcessed) ;
708 //PostData(2,fCFManager->GetParticleContainer()) ;
709 //PostData(3,fCorrelation) ;
713 //___________________________________________________________________________
714 void AliCFHeavyFlavourTaskMultiVarMultiStep::Terminate(Option_t*)
716 // The Terminate() function is the last function to be called during
717 // a query. It always runs on the client, it can be used to present
718 // the results graphically or save the results to file.
720 AliAnalysisTaskSE::Terminate();
722 AliInfo(Form("Found %i MC particles that are D0 in MC, in %d events",fCountMC,fEvents));
723 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, in %d events",fCountAcc,fEvents));
724 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fEvents));
725 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));
726 AliInfo(Form("Found %i reco D0 that are decaying in K+pi, in %d events",fCountReco,fEvents));
727 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));
728 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));
729 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and satisfy PPR cuts, in %d events",fCountRecoPPR,fEvents));
731 // draw some example plots....
733 // AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
734 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
736 printf("CONTAINER NOT FOUND\n");
739 // projecting the containers to obtain histograms
740 // first argument = variable, second argument = step
743 TH1D* h00 = cont->ShowProjection(0,0) ; // pt
744 TH1D* h10 = cont->ShowProjection(1,0) ; // rapidity
745 TH1D* h20 = cont->ShowProjection(2,0) ; // cosThetaStar
746 TH1D* h30 = cont->ShowProjection(3,0) ; // pTpi
747 TH1D* h40 = cont->ShowProjection(4,0) ; // pTK
748 TH1D* h50 = cont->ShowProjection(5,0) ; // cT
749 TH1D* h60 = cont->ShowProjection(6,0) ; // dca
750 TH1D* h70 = cont->ShowProjection(7,0) ; // d0pi
751 TH1D* h80 = cont->ShowProjection(8,0) ; // d0K
752 TH1D* h90 = cont->ShowProjection(9,0) ; // d0xd0
753 TH1D* h100 = cont->ShowProjection(10,0) ; // cosPointingAngle
754 TH1D* h110 = cont->ShowProjection(11,0) ; // phi
756 // MC-Acceptance level
757 TH1D* h01 = cont->ShowProjection(0,1) ; // pt
758 TH1D* h11 = cont->ShowProjection(1,1) ; // rapidity
759 TH1D* h21 = cont->ShowProjection(2,1) ; // cosThetaStar
760 TH1D* h31 = cont->ShowProjection(3,1) ; // pTpi
761 TH1D* h41 = cont->ShowProjection(4,1) ; // pTK
762 TH1D* h51 = cont->ShowProjection(5,1) ; // cT
763 TH1D* h61 = cont->ShowProjection(6,1) ; // dca
764 TH1D* h71 = cont->ShowProjection(7,1) ; // d0pi
765 TH1D* h81 = cont->ShowProjection(8,1) ; // d0K
766 TH1D* h91 = cont->ShowProjection(9,1) ; // d0xd0
767 TH1D* h101 = cont->ShowProjection(10,1) ; // cosPointingAngle
768 TH1D* h111 = cont->ShowProjection(11,1) ; // phi
771 TH1D* h04 = cont->ShowProjection(0,4) ; // pt
772 TH1D* h14 = cont->ShowProjection(1,4) ; // rapidity
773 TH1D* h24 = cont->ShowProjection(2,4) ; // cosThetaStar
774 TH1D* h34 = cont->ShowProjection(3,4) ; // pTpi
775 TH1D* h44 = cont->ShowProjection(4,4) ; // pTK
776 TH1D* h54 = cont->ShowProjection(5,4) ; // cT
777 TH1D* h64 = cont->ShowProjection(6,4) ; // dca
778 TH1D* h74 = cont->ShowProjection(7,4) ; // d0pi
779 TH1D* h84 = cont->ShowProjection(8,4) ; // d0K
780 TH1D* h94 = cont->ShowProjection(9,4) ; // d0xd0
781 TH1D* h104 = cont->ShowProjection(10,4) ; // cosPointingAngle
782 TH1D* h114 = cont->ShowProjection(11,4) ; // phi
784 h00->SetTitle("pT_D0 (GeV/c)");
785 h10->SetTitle("rapidity");
786 h20->SetTitle("cosThetaStar");
787 h30->SetTitle("pT_pi (GeV/c)");
788 h40->SetTitle("pT_K (Gev/c)");
789 h50->SetTitle("cT (#mum)");
790 h60->SetTitle("dca (#mum)");
791 h70->SetTitle("d0_pi (#mum)");
792 h80->SetTitle("d0_K (#mum)");
793 h90->SetTitle("d0xd0 (#mum^2)");
794 h100->SetTitle("cosPointingAngle");
795 h100->SetTitle("phi (rad)");
797 h00->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
798 h10->GetXaxis()->SetTitle("rapidity");
799 h20->GetXaxis()->SetTitle("cosThetaStar");
800 h30->GetXaxis()->SetTitle("pT_pi (GeV/c)");
801 h40->GetXaxis()->SetTitle("pT_K (Gev/c)");
802 h50->GetXaxis()->SetTitle("cT (#mum)");
803 h60->GetXaxis()->SetTitle("dca (#mum)");
804 h70->GetXaxis()->SetTitle("d0_pi (#mum)");
805 h80->GetXaxis()->SetTitle("d0_K (#mum)");
806 h90->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
807 h100->GetXaxis()->SetTitle("cosPointingAngle");
808 h110->GetXaxis()->SetTitle("phi (rad)");
810 h01->SetTitle("pT_D0 (GeV/c)");
811 h11->SetTitle("rapidity");
812 h21->SetTitle("cosThetaStar");
813 h31->SetTitle("pT_pi (GeV/c)");
814 h41->SetTitle("pT_K (Gev/c)");
815 h51->SetTitle("cT (#mum)");
816 h61->SetTitle("dca (#mum)");
817 h71->SetTitle("d0_pi (#mum)");
818 h81->SetTitle("d0_K (#mum)");
819 h91->SetTitle("d0xd0 (#mum^2)");
820 h101->SetTitle("cosPointingAngle");
821 h111->GetXaxis()->SetTitle("phi (rad)");
823 h01->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
824 h11->GetXaxis()->SetTitle("rapidity");
825 h21->GetXaxis()->SetTitle("cosThetaStar");
826 h31->GetXaxis()->SetTitle("pT_pi (GeV/c)");
827 h41->GetXaxis()->SetTitle("pT_K (Gev/c)");
828 h51->GetXaxis()->SetTitle("cT (#mum)");
829 h61->GetXaxis()->SetTitle("dca (#mum)");
830 h71->GetXaxis()->SetTitle("d0_pi (#mum)");
831 h81->GetXaxis()->SetTitle("d0_K (#mum)");
832 h91->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
833 h101->GetXaxis()->SetTitle("cosPointingAngle");
834 h111->GetXaxis()->SetTitle("phi (rad)");
836 h04->SetTitle("pT_D0 (GeV/c)");
837 h14->SetTitle("rapidity");
838 h24->SetTitle("cosThetaStar");
839 h34->SetTitle("pT_pi (GeV/c)");
840 h44->SetTitle("pT_K (Gev/c)");
841 h54->SetTitle("cT (#mum)");
842 h64->SetTitle("dca (#mum)");
843 h74->SetTitle("d0_pi (#mum)");
844 h84->SetTitle("d0_K (#mum)");
845 h94->SetTitle("d0xd0 (#mum^2)");
846 h104->SetTitle("cosPointingAngle");
847 h114->GetXaxis()->SetTitle("phi (rad)");
849 h04->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
850 h14->GetXaxis()->SetTitle("rapidity");
851 h24->GetXaxis()->SetTitle("cosThetaStar");
852 h34->GetXaxis()->SetTitle("pT_pi (GeV/c)");
853 h44->GetXaxis()->SetTitle("pT_K (Gev/c)");
854 h54->GetXaxis()->SetTitle("cT (#mum)");
855 h64->GetXaxis()->SetTitle("dca (#mum)");
856 h74->GetXaxis()->SetTitle("d0_pi (#mum)");
857 h84->GetXaxis()->SetTitle("d0_K (#mum)");
858 h94->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
859 h104->GetXaxis()->SetTitle("cosPointingAngle");
860 h114->GetXaxis()->SetTitle("phi (rad)");
862 Double_t max0 = h00->GetMaximum();
863 Double_t max1 = h10->GetMaximum();
864 Double_t max2 = h20->GetMaximum();
865 Double_t max3 = h30->GetMaximum();
866 Double_t max4 = h40->GetMaximum();
867 Double_t max5 = h50->GetMaximum();
868 Double_t max6 = h60->GetMaximum();
869 Double_t max7 = h70->GetMaximum();
870 Double_t max8 = h80->GetMaximum();
871 Double_t max9 = h90->GetMaximum();
872 Double_t max10 = h100->GetMaximum();
873 Double_t max11 = h110->GetMaximum();
875 h00->GetYaxis()->SetRangeUser(0,max0*1.2);
876 h10->GetYaxis()->SetRangeUser(0,max1*1.2);
877 h20->GetYaxis()->SetRangeUser(0,max2*1.2);
878 h30->GetYaxis()->SetRangeUser(0,max3*1.2);
879 h40->GetYaxis()->SetRangeUser(0,max4*1.2);
880 h50->GetYaxis()->SetRangeUser(0,max5*1.2);
881 h60->GetYaxis()->SetRangeUser(0,max6*1.2);
882 h70->GetYaxis()->SetRangeUser(0,max7*1.2);
883 h80->GetYaxis()->SetRangeUser(0,max8*1.2);
884 h90->GetYaxis()->SetRangeUser(0,max9*1.2);
885 h100->GetYaxis()->SetRangeUser(0,max10*1.2);
886 h110->GetYaxis()->SetRangeUser(0,max11*1.2);
888 h01->GetYaxis()->SetRangeUser(0,max0*1.2);
889 h11->GetYaxis()->SetRangeUser(0,max1*1.2);
890 h21->GetYaxis()->SetRangeUser(0,max2*1.2);
891 h31->GetYaxis()->SetRangeUser(0,max3*1.2);
892 h41->GetYaxis()->SetRangeUser(0,max4*1.2);
893 h51->GetYaxis()->SetRangeUser(0,max5*1.2);
894 h61->GetYaxis()->SetRangeUser(0,max6*1.2);
895 h71->GetYaxis()->SetRangeUser(0,max7*1.2);
896 h81->GetYaxis()->SetRangeUser(0,max8*1.2);
897 h91->GetYaxis()->SetRangeUser(0,max9*1.2);
898 h101->GetYaxis()->SetRangeUser(0,max10*1.2);
899 h111->GetYaxis()->SetRangeUser(0,max11*1.2);
901 h00->SetMarkerStyle(20);
902 h10->SetMarkerStyle(24);
903 h20->SetMarkerStyle(21);
904 h30->SetMarkerStyle(25);
905 h40->SetMarkerStyle(27);
906 h50->SetMarkerStyle(28);
907 h60->SetMarkerStyle(20);
908 h70->SetMarkerStyle(24);
909 h80->SetMarkerStyle(21);
910 h90->SetMarkerStyle(25);
911 h100->SetMarkerStyle(27);
912 h110->SetMarkerStyle(28);
914 h00->SetMarkerColor(2);
915 h10->SetMarkerColor(2);
916 h20->SetMarkerColor(2);
917 h30->SetMarkerColor(2);
918 h40->SetMarkerColor(2);
919 h50->SetMarkerColor(2);
920 h60->SetMarkerColor(2);
921 h70->SetMarkerColor(2);
922 h80->SetMarkerColor(2);
923 h90->SetMarkerColor(2);
924 h100->SetMarkerColor(2);
925 h110->SetMarkerColor(2);
927 h01->SetMarkerStyle(20) ;
928 h11->SetMarkerStyle(24) ;
929 h21->SetMarkerStyle(21) ;
930 h31->SetMarkerStyle(25) ;
931 h41->SetMarkerStyle(27) ;
932 h51->SetMarkerStyle(28) ;
933 h61->SetMarkerStyle(20);
934 h71->SetMarkerStyle(24);
935 h81->SetMarkerStyle(21);
936 h91->SetMarkerStyle(25);
937 h101->SetMarkerStyle(27);
938 h111->SetMarkerStyle(28);
940 h01->SetMarkerColor(8);
941 h11->SetMarkerColor(8);
942 h21->SetMarkerColor(8);
943 h31->SetMarkerColor(8);
944 h41->SetMarkerColor(8);
945 h51->SetMarkerColor(8);
946 h61->SetMarkerColor(8);
947 h71->SetMarkerColor(8);
948 h81->SetMarkerColor(8);
949 h91->SetMarkerColor(8);
950 h101->SetMarkerColor(8);
951 h111->SetMarkerColor(8);
953 h04->SetMarkerStyle(20) ;
954 h14->SetMarkerStyle(24) ;
955 h24->SetMarkerStyle(21) ;
956 h34->SetMarkerStyle(25) ;
957 h44->SetMarkerStyle(27) ;
958 h54->SetMarkerStyle(28) ;
959 h64->SetMarkerStyle(20);
960 h74->SetMarkerStyle(24);
961 h84->SetMarkerStyle(21);
962 h94->SetMarkerStyle(25);
963 h104->SetMarkerStyle(27);
964 h114->SetMarkerStyle(28);
966 h04->SetMarkerColor(4);
967 h14->SetMarkerColor(4);
968 h24->SetMarkerColor(4);
969 h34->SetMarkerColor(4);
970 h44->SetMarkerColor(4);
971 h54->SetMarkerColor(4);
972 h64->SetMarkerColor(4);
973 h74->SetMarkerColor(4);
974 h84->SetMarkerColor(4);
975 h94->SetMarkerColor(4);
976 h104->SetMarkerColor(4);
977 h114->SetMarkerColor(4);
979 gStyle->SetCanvasColor(0);
980 gStyle->SetFrameFillColor(0);
981 gStyle->SetTitleFillColor(0);
982 gStyle->SetStatColor(0);
984 // drawing in 2 separate canvas for a matter of clearity
985 TCanvas * c1 =new TCanvas("c1","pT, rapidiy, cosThetaStar",1100,1600);
1017 TCanvas * c2 =new TCanvas("c2","pTpi, pTK, cT",1100,1600);
1048 TCanvas * c3 =new TCanvas("c3","dca, d0pi, d0K",1100,1600);
1079 TCanvas * c4 =new TCanvas("c4","d0xd0, cosPointingAngle, phi",1100,1600);
1110 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
1112 TH2D* corr1 =hcorr->Projection(0,2);
1113 TH2D* corr2 = hcorr->Projection(1,3);
1115 TCanvas * c7 =new TCanvas("c7","",800,400);
1118 corr1->Draw("text");
1120 corr2->Draw("text");
1123 TString projection_filename="CFtaskHFprojection";
1124 if(fKeepD0fromBOnly) projection_filename+="_KeepD0fromBOnly";
1125 else if(fKeepD0fromB) projection_filename+="_KeepD0fromB";
1126 projection_filename+=".root";
1127 TFile* file_projection = new TFile(projection_filename.Data(),"RECREATE");
1131 h00->Write("pT_D0_step0");
1132 h10->Write("rapidity_step0");
1133 h20->Write("cosThetaStar_step0");
1134 h30->Write("pT_pi_step0");
1135 h40->Write("pT_K_step0");
1136 h50->Write("cT_step0");
1137 h60->Write("dca_step0");
1138 h70->Write("d0_pi_step0");
1139 h80->Write("d0_K_step0");
1140 h90->Write("d0xd0_step0");
1141 h100->Write("cosPointingAngle_step0");
1142 h110->Write("phi_step0");
1144 h01->Write("pT_D0_step1");
1145 h11->Write("rapidity_step1");
1146 h21->Write("cosThetaStar_step1");
1147 h31->Write("pT_pi_step1");
1148 h41->Write("pT_K_step1");
1149 h51->Write("cT_step1");
1150 h61->Write("dca_step1");
1151 h71->Write("d0_pi_step1");
1152 h81->Write("d0_K_step1");
1153 h91->Write("d0xd0_step1");
1154 h101->Write("cosPointingAngle_step1");
1155 h111->Write("phi_step1");
1157 h04->Write("pT_D0_step2");
1158 h14->Write("rapidity_step2");
1159 h24->Write("cosThetaStar_step2");
1160 h34->Write("pT_pi_step2");
1161 h44->Write("pT_K_step2");
1162 h54->Write("cT_step2");
1163 h64->Write("dca_step2");
1164 h74->Write("d0_pi_step2");
1165 h80->Write("d0_K_step2");
1166 h94->Write("d0xd0_step2");
1167 h104->Write("cosPointingAngle_step2");
1168 h114->Write("phi_step2");
1170 file_projection->Close();
1175 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
1176 c2->SaveAs("Plots/pTpi_pTK_cT.eps");
1177 c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
1178 c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
1180 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
1181 c2->SaveAs("Plots/pTpi_pTK_cT.gif");
1182 c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
1183 c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
1187 //___________________________________________________________________________
1188 void AliCFHeavyFlavourTaskMultiVarMultiStep::UserCreateOutputObjects() {
1189 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1190 //TO BE SET BEFORE THE EXECUTION OF THE TASK
1192 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1196 fHistEventsProcessed = new TH1I("CFHFchist0","",1,0,1) ;
1199 //___________________________________________________________________________
1200 Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CosThetaStar(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
1203 // to calculate cos(ThetaStar) for generated particle
1204 // using the K, since mcPartDaughter0 and mcPartDaughter1 always correspond to K and pi respectively
1205 // (see where the function is called)
1208 Int_t pdgvtx = mcPart->GetPdgCode();
1209 /* if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1210 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1211 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1212 AliInfo(Form("D0, with pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1213 AliDebug(2,"This is a D0");
1214 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1215 mcPartDaughter0 = mcPartDaughter1;
1216 mcPartDaughter1 = mcPartdummy;
1222 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1223 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1224 if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1228 AliDebug(2,"D0bar");
1230 if (pdgprong0 == 211){
1231 AliDebug(2,Form("pdgprong0 = %d, pdgprong1 = %d, switching...",pdgprong0,pdgprong1));
1232 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1233 mcPartDaughter0 = mcPartDaughter1;
1234 mcPartDaughter1 = mcPartdummy;
1235 pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1236 pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1239 AliDebug(2,Form("After checking, pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1240 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdgvtx))->Mass();
1242 massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
1243 massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
1245 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);
1247 Double_t px = mcPartDaughter0->Px()+mcPartDaughter1->Px();
1248 Double_t py = mcPartDaughter0->Py()+mcPartDaughter1->Py();
1249 Double_t pz = mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
1250 Double_t p = TMath::Sqrt(px*px+py*py+pz*pz);
1251 Double_t e = TMath::Sqrt(massvtx*massvtx+p*p);
1252 Double_t beta = p/e;
1253 Double_t gamma = e/massvtx;
1254 TVector3 mom(mcPartDaughter0->Px(),mcPartDaughter0->Py(),mcPartDaughter0->Pz());
1255 TVector3 momTot(mcPartDaughter0->Px()+mcPartDaughter1->Px(),mcPartDaughter0->Py()+mcPartDaughter1->Py(),mcPartDaughter0->Pz()+mcPartDaughter1->Pz());
1257 Double_t qlprong = mom.Dot(momTot)/momTot.Mag(); // analog to AliAODRecoDecay::QlProng(0)
1259 AliDebug(2,Form("pStar = %f, beta = %f, gamma = %f, qlprong = %f, massp[0] = %f", pStar, beta, gamma, qlprong, massp[0]));
1260 Double_t cts = (qlprong/gamma-beta*TMath::Sqrt(pStar*pStar+massp[0]*massp[0]))/pStar;
1261 AliDebug(2,Form("cts = %f", cts));
1264 //___________________________________________________________________________
1265 Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CT(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
1268 // to calculate cT for generated particle
1271 Double_t xmcPart[3] = {0,0,0};
1272 Double_t xdaughter0[3] = {0,0,0};
1273 Double_t xdaughter1[3] = {0,0,0};
1274 mcPart->XvYvZv(xmcPart); // cm
1275 mcPartDaughter0->XvYvZv(xdaughter0); // cm
1276 mcPartDaughter1->XvYvZv(xdaughter1); //cm
1277 Double_t prodVtxD0 = TMath::Sqrt(xmcPart[0]*xmcPart[0]+
1278 xmcPart[1]*xmcPart[1]+
1279 xmcPart[2]*xmcPart[2]);
1280 Double_t prodVtxDaughter0 = TMath::Sqrt(xdaughter0[0]*xdaughter0[0]+
1281 xdaughter0[1]*xdaughter0[1]+
1282 xdaughter0[2]*xdaughter0[2]);
1283 Double_t prodVtxDaughter1 = TMath::Sqrt(xdaughter1[0]*xdaughter1[0]+
1284 xdaughter1[1]*xdaughter1[1]+
1285 xdaughter1[2]*xdaughter1[2]);
1287 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));
1288 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));
1289 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));
1291 Double_t cT0 = TMath::Sqrt((xdaughter0[0]-xmcPart[0])*(xdaughter0[0]-xmcPart[0])+
1292 (xdaughter0[1]-xmcPart[1])*(xdaughter0[1]-xmcPart[1])+
1293 (xdaughter0[2]-xmcPart[2])*(xdaughter0[2]-xmcPart[2]));
1295 Double_t cT1 = TMath::Sqrt((xdaughter1[0]-xmcPart[0])*(xdaughter1[0]-xmcPart[0])+
1296 (xdaughter1[1]-xmcPart[1])*(xdaughter1[1]-xmcPart[1])+
1297 (xdaughter1[2]-xmcPart[2])*(xdaughter1[2]-xmcPart[2]));
1299 if ((cT0 - cT1)>1E-5) {
1300 AliWarning(Form("cT from daughter 0 (%f) different from cT from daughter 1 (%f)! Using cT from daughter 0, but PLEASE, CHECK....",cT0,cT1));
1303 // calculating cT from cT0
1305 Double_t mass = TDatabasePDG::Instance()->GetParticle(mcPart->GetPdgCode())->Mass(); // mcPart->GetPdgCode() should return +/- 421 for the D0/D0bar
1306 Double_t p = mcPart-> P();
1307 Double_t cT = cT0*mass/p;
1308 AliDebug(2, Form("cT from daughter 0 = %f (micron)", cT0*1E4));
1309 AliDebug(2, Form("cT from daughter 1 = %f (micron)", cT1*1E4));
1310 AliDebug(2, Form("cT (with mass = %f and p = %f) = %f (micron)", mass, p, cT*1E4));
1313 //________________________________________________________________________________
1314 Bool_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC) const {
1317 // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
1320 Bool_t isOk = kFALSE;
1322 // check whether the D0 decays in pi+K
1323 // to be added!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1324 // could use a cut in the CF, but not clear how to define a TDecayChannel
1325 // implemented for the time being as a cut on the number of daughters - see later when
1326 // getting the daughters
1328 // getting the daughters
1329 Int_t daughter0 = mcPart->GetDaughter(0);
1330 Int_t daughter1 = mcPart->GetDaughter(1);
1331 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
1332 if (daughter0 == 0 || daughter1 == 0) {
1333 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
1336 if (TMath::Abs(daughter1 - daughter0) != 1) {
1337 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
1340 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
1341 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
1342 if (!mcPartDaughter0 || !mcPartDaughter1) {
1343 AliWarning("At least one Daughter Particle not found in tree, skipping");
1346 if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==321 &&
1347 TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
1348 !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
1349 TMath::Abs(mcPartDaughter1->GetPdgCode())==321)) {
1350 AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
1354 Double_t vtx1[3] = {0,0,0}; // primary vertex
1355 Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
1356 Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
1357 mcPart->XvYvZv(vtx1); // cm
1358 // getting vertex from daughters
1359 mcPartDaughter0->XvYvZv(vtx2daughter0); // cm
1360 mcPartDaughter1->XvYvZv(vtx2daughter1); //cm
1361 if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
1362 AliError("Daughters have different secondary vertex, skipping the track");
1367 // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
1368 AliAODMCParticle* positiveDaugh = mcPartDaughter0;
1369 AliAODMCParticle* negativeDaugh = mcPartDaughter1;
1370 if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
1371 // inverting in case the positive daughter is the second one
1372 positiveDaugh = mcPartDaughter1;
1373 negativeDaugh = mcPartDaughter0;
1375 // getting the momentum from the daughters
1376 Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
1377 Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
1378 Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
1380 Double_t d0[2] = {0.,0.};
1382 AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
1384 Double_t cosThetaStar = 0.;
1385 Double_t cosThetaStarD0 = 0.;
1386 Double_t cosThetaStarD0bar = 0.;
1387 cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
1388 cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
1389 if (mcPart->GetPdgCode() == 421){ // D0
1390 AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1391 cosThetaStar = cosThetaStarD0;
1393 else if (mcPart->GetPdgCode() == -421){ // D0bar{
1394 AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1395 cosThetaStar = cosThetaStarD0bar;
1398 AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
1401 if (cosThetaStar < -1 || cosThetaStar > 1) {
1402 AliWarning("Invalid value for cosine Theta star, returning");
1406 // calculate cos(Theta*) according to the method implemented herein
1408 Double_t cts = 9999.;
1409 cts = CosThetaStar(mcPart, mcPartDaughter0, mcPartDaughter1);
1410 if (cts < -1 || cts > 1) {
1411 AliWarning("Invalid value for cosine Theta star from AliCFHeavyFlavourTaskMultiVarMultiStep method");
1413 if (TMath::Abs(cts - cosThetaStar)>0.001) {
1414 AliError(Form("cosThetaStar automatically calculated different from that manually calculated!!! cosThetaStar = %f, cosThetaStar = %f", cosThetaStar,cts));
1417 Double_t cT = decay->Ct(421);
1419 // calculate cT from the method implemented herein
1421 cT1 = CT(mcPart, mcPartDaughter0, mcPartDaughter1);
1423 if (TMath::Abs(cT1 - cT)>0.001) {
1424 AliError(Form("cT automatically calculated different from that manually calculated!!! cT = %f, cT1 = %f",cT,cT1));
1427 // get the pT of the daughters
1432 if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
1433 pTpi = mcPartDaughter0->Pt();
1434 pTK = mcPartDaughter1->Pt();
1437 pTpi = mcPartDaughter1->Pt();
1438 pTK = mcPartDaughter0->Pt();
1441 vectorMC[0] = mcPart->Pt();
1442 vectorMC[1] = mcPart->Y() ;
1443 vectorMC[2] = cosThetaStar ;
1444 vectorMC[3] = pTpi ;
1446 vectorMC[5] = cT*1.E4 ; // in micron
1447 vectorMC[6] = mcPart->Phi() ;
1451 //_________________________________________________________________________________________________
1452 Int_t AliCFHeavyFlavourTaskMultiVarMultiStep::CheckOrigin(AliAODMCParticle* mcPart, TClonesArray* mcArray)const{
1455 // checking whether the very mother of the D0 is a charm or a bottom quark
1458 Int_t pdgGranma = 0;
1460 mother = mcPart->GetMother();
1464 AliDebug(2,Form("mother at step %d = %d", istep, mother));
1465 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
1466 pdgGranma = mcGranma->GetPdgCode();
1467 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1468 Int_t abspdgGranma = TMath::Abs(pdgGranma);
1469 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1472 mother = mcGranma->GetMother();
1476 //__________________________________________________________________________________________________
1477 Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetWeight(Float_t pt){
1480 // calculating the weight to fill the container
1507 Double_t func1[4] = {0.2149,2.433,2.102,2.5};
1508 Double_t func2[4] = {3.63832E-2,1.88238,1.34854,2.5};
1510 Double_t dndpt_func1 = dNdptFit(pt,func1);
1511 Double_t dndpt_func2 = dNdptFit(pt,func2);
1512 AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
1513 return dndpt_func1/dndpt_func2;
1516 //__________________________________________________________________________________________________
1517 Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::dNdptFit(Float_t pt, Double_t* par){
1520 // calculating dNdpt
1523 Double_t denom = TMath::Power((pt/par[1]), par[3] );
1524 Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);