1 /**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 //-----------------------------------------------------------------------
17 // Class for HF corrections as a function of many variables
18 // 6 Steps introduced: MC, MC Acc, Reco, Reco Acc, Reco Acc + ITS Cl,
19 // Reco Acc + ITS Cl + PPR cuts
20 // 12 variables used: pt, y, cosThetaStar, ptPi, ptK, ct,
21 // dca, d0Pi, d0K, d0Pixd0K, cosPointingAngle, phi
23 //-----------------------------------------------------------------------
24 // Author : C. Zampolli, CERN
25 //-----------------------------------------------------------------------
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 "AliAODEvent.h"
45 #include "AliAODRecoDecay.h"
46 #include "AliAODRecoDecayHF.h"
47 #include "AliAODRecoDecayHF2Prong.h"
48 #include "AliAODMCParticle.h"
49 #include "AliESDtrack.h"
51 #include "THnSparse.h"
53 //__________________________________________________________________________
54 AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep() :
58 fHistEventsProcessed(0x0),
66 fCountRecoITSClusters(0),
69 fFillFromGenerated(kFALSE),
77 //___________________________________________________________________________
78 AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(const Char_t* name) :
79 AliAnalysisTaskSE(name),
82 fHistEventsProcessed(0x0),
90 fCountRecoITSClusters(0),
93 fFillFromGenerated(kFALSE),
98 // Constructor. Initialization of Inputs and Outputs
100 Info("AliCFHeavyFlavourTaskMultiVarMultiStep","Calling Constructor");
102 DefineInput(0) and DefineOutput(0)
103 are taken care of by AliAnalysisTaskSE constructor
105 DefineOutput(1,TH1I::Class());
106 DefineOutput(2,AliCFContainer::Class());
107 DefineOutput(3,THnSparseD::Class());
110 //___________________________________________________________________________
111 AliCFHeavyFlavourTaskMultiVarMultiStep& AliCFHeavyFlavourTaskMultiVarMultiStep::operator=(const AliCFHeavyFlavourTaskMultiVarMultiStep& c)
114 // Assignment operator
117 AliAnalysisTaskSE::operator=(c) ;
119 fCFManager = c.fCFManager;
120 fHistEventsProcessed = c.fHistEventsProcessed;
125 //___________________________________________________________________________
126 AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(const AliCFHeavyFlavourTaskMultiVarMultiStep& c) :
127 AliAnalysisTaskSE(c),
129 fCFManager(c.fCFManager),
130 fHistEventsProcessed(c.fHistEventsProcessed),
131 fCorrelation(c.fCorrelation),
132 fCountMC(c.fCountMC),
133 fCountAcc(c.fCountAcc),
134 fCountVertex(c.fCountVertex),
135 fCountRefit(c.fCountRefit),
136 fCountReco(c.fCountReco),
137 fCountRecoAcc(c.fCountRecoAcc),
138 fCountRecoITSClusters(c.fCountRecoITSClusters),
139 fCountRecoPPR(c.fCountRecoPPR),
141 fFillFromGenerated(c.fFillFromGenerated),
142 fMinITSClusters(c.fMinITSClusters),
143 fAcceptanceUnf(c.fAcceptanceUnf)
150 //___________________________________________________________________________
151 AliCFHeavyFlavourTaskMultiVarMultiStep::~AliCFHeavyFlavourTaskMultiVarMultiStep() {
155 if (fCFManager) delete fCFManager ;
156 if (fHistEventsProcessed) delete fHistEventsProcessed ;
157 if (fCorrelation) delete fCorrelation ;
160 //_________________________________________________
161 void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
164 // Main loop function
167 if (fFillFromGenerated){
168 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
172 Error("UserExec","NO EVENT FOUND!");
177 if (fEvents%10000 ==0) AliDebug(2,Form("Event %d",fEvents));
178 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
179 fCFManager->SetEventInfo(aodEvent);
181 // MC-event selection
182 Double_t containerInput[12] ;
183 Double_t containerInputMC[12] ;
185 //loop on the MC event
187 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
188 if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
191 Int_t icountReco = 0;
192 Int_t icountVertex = 0;
193 Int_t icountRefit = 0;
194 Int_t icountRecoAcc = 0;
195 Int_t icountRecoITSClusters = 0;
196 Int_t icountRecoPPR = 0;
200 // AOD primary vertex
201 AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
202 Bool_t vtxFlag = kTRUE;
203 TString title=vtx1->GetTitle();
204 if(!title.Contains("VertexerTracks")) vtxFlag=kFALSE;
206 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
207 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
208 if (mcPart->GetPdgCode() == 4) cquarks++;
209 if (mcPart->GetPdgCode() == -4) cquarks++;
211 AliWarning("Particle not found in tree, skipping");
215 // check the MC-level cuts
217 if (!fCFManager->CheckParticleCuts(0, mcPart)) continue; // 0 stands for MC level
218 Int_t pdgGranma = CheckOrigin(mcPart, mcArray);
219 Int_t abspdgGranma = TMath::Abs(pdgGranma);
220 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
221 AliDebug(2,Form("Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
222 continue; // skipping particles that don't come from c quark
225 // if (TMath::Abs(pdgGranma)!=4) {
227 // fill the container for Gen-level selection
228 Double_t vectorMC[7] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.};
229 if (GetGeneratedValuesFromMCParticle(mcPart, mcArray, vectorMC)){
230 containerInputMC[0] = vectorMC[0];
231 containerInputMC[1] = vectorMC[1] ;
232 containerInputMC[2] = vectorMC[2] ;
233 containerInputMC[3] = vectorMC[3] ;
234 containerInputMC[4] = vectorMC[4] ;
235 containerInputMC[5] = vectorMC[5] ; // in micron
236 containerInputMC[6] = 0.; // dummy value, meaningless in MC, in micron
237 containerInputMC[7] = 0.; // dummy value, meaningless in MC, in micron
238 containerInputMC[8] = 0.; // dummy value, meaningless in MC, in micron
239 containerInputMC[9] = -100000.; // dummy value, meaningless in MC, in micron^2
240 containerInputMC[10] = 1.01; // dummy value, meaningless in MC
241 containerInputMC[11] = vectorMC[6]; // dummy value, meaningless in MC
242 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGenerated);
245 // check the MC-Acceptance level cuts
246 // since standard CF functions are not applicable, using Kine Cuts on daughters
248 Int_t daughter0 = mcPart->GetDaughter(0);
249 Int_t daughter1 = mcPart->GetDaughter(1);
250 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
251 if (daughter0 == 0 || daughter1 == 0) {
252 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!! But it should have, this check was already done...");
254 if (TMath::Abs(daughter1 - daughter0) != 1) {
255 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, but it should be, this check was already done...");
257 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
258 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
259 if (!mcPartDaughter0 || !mcPartDaughter1) {
260 AliWarning("At least one Daughter Particle not found in tree, but it should be, this check was already done...");
262 Double_t eta0 = mcPartDaughter0->Eta();
263 Double_t eta1 = mcPartDaughter1->Eta();
264 Double_t y0 = mcPartDaughter0->Y();
265 Double_t y1 = mcPartDaughter1->Y();
266 Double_t pt0 = mcPartDaughter0->Pt();
267 Double_t pt1 = mcPartDaughter1->Pt();
268 AliDebug(2, Form("Daughter 0: eta = %f, y = %f, pt = %f", eta0, y0, pt0));
269 AliDebug(2, Form("Daughter 1: eta = %f, y = %f, pt = %f", eta1, y1, pt1));
270 Bool_t daught0inAcceptance = (TMath::Abs(eta0) < 0.9 && pt0 > 0.1);
271 Bool_t daught1inAcceptance = (TMath::Abs(eta1) < 0.9 && pt1 > 0.1);
272 if (daught0inAcceptance && daught1inAcceptance) {
273 // checking whether the cuts implemented in the CF are equivalent - simply a cross-check
274 AliDebug(2, "Daughter particles in acceptance");
275 if (!fCFManager->CheckParticleCuts(1, mcPartDaughter0)) {
276 AliDebug(2,"Inconsistency with CF cut for daughter 0!");
278 if (!fCFManager->CheckParticleCuts(1, mcPartDaughter1)) {
279 AliDebug(2,"Inconsistency with CF cut for daughter 1!");
281 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance);
284 // check on the vertex
286 printf("Vertex cut passed\n");
287 // filling the container if the vertex is ok
288 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex) ;
289 printf("Vertex cut passed 2\n");
291 // check on the kTPCrefit and kITSrefit conditions of the daughters
292 Bool_t refitFlag = kTRUE;
293 for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
294 AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
295 if ((track->GetLabel() == daughter0) || (track->GetLabel() == daughter1)) {
296 if(!(track->GetStatus()&AliESDtrack::kTPCrefit) || !(track->GetStatus()&AliESDtrack::kITSrefit)) {
302 printf("Refit cut passed\n");
303 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit);
307 AliDebug(3,"Refit cut not passed\n");
311 AliDebug(3,"Vertex cut not passed\n");
315 AliDebug(3,"Acceptance cut not passed\n");
319 AliDebug(3,"Problems in filling the container");
324 if (cquarks<2) AliDebug(2,Form("Event found with %d c-quarks", cquarks));
326 AliDebug(2,Form("Found %i MC particles that are D0!!",icountMC));
327 AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Acc cuts!!",icountAcc));
328 AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Vertex cuts!!",icountVertex));
329 AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Refit cuts!!",icountRefit));
331 // Now go to rec level
332 fCountMC += icountMC;
333 fCountAcc += icountAcc;
335 // load heavy flavour vertices
336 TClonesArray *arrayD0toKpi = (TClonesArray*)((aodEvent->GetList())->FindObject("D0toKpi"));
337 if (!arrayD0toKpi) AliError("Could not find array of HF vertices");
338 AliDebug(2, Form("Found %d vertices",arrayD0toKpi->GetEntriesFast()));
340 Int_t pdgDgD0toKpi[2]={321,211};
342 for (Int_t iD0toKpi = 0; iD0toKpi<arrayD0toKpi->GetEntriesFast(); iD0toKpi++) {
344 AliAODRecoDecayHF2Prong* d0tokpi = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->At(iD0toKpi);
345 Bool_t unsetvtx=kFALSE;
346 if(!d0tokpi->GetOwnPrimaryVtx()) {
347 d0tokpi->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
351 // find associated MC particle
352 Int_t mcLabel = d0tokpi->MatchToMC(421,mcArray,2,pdgDgD0toKpi) ;
355 AliDebug(2,"No MC particle found");
358 AliAODMCParticle* mcVtxHF = (AliAODMCParticle*)mcArray->At(mcLabel);
360 AliWarning("Could not find associated MC in AOD MC tree");
363 // check whether the daughters have kTPCrefit and kITSrefit set
364 AliAODTrack *track0 = (AliAODTrack*)d0tokpi->GetDaughter(0);
365 AliAODTrack *track1 = (AliAODTrack*)d0tokpi->GetDaughter(1);
366 if((!(track0->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track1->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track0->GetStatus()&AliESDtrack::kITSrefit)) || (!(track1->GetStatus()&AliESDtrack::kITSrefit))) {
367 // skipping if at least one daughter does not have kTPCrefit or kITSrefit
371 // check on the vertex -- could be moved outside the loop on the reconstructed D0...
373 // skipping cause vertex is not ok
377 const Double_t d0tokpiCuts[9] = {0.3,999999.,1.1,0.,0.,999999.,999999.,999999.,0.};
379 if (!(d0tokpi->SelectD0(&d0tokpiCuts[0],okD0,okD0bar))){
380 // skipping candidate
384 // check if associated MC v0 passes the cuts
385 if (!fCFManager->CheckParticleCuts(0 ,mcVtxHF)) { // 0 stands for MC
386 AliDebug(2, "Skipping the particles due to cuts");
389 Int_t pdgGranma = CheckOrigin(mcVtxHF, mcArray);
390 Int_t abspdgGranma = TMath::Abs(pdgGranma);
391 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
392 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));
393 continue; // skipping particles that don't come from c quark
396 // fill the container...
397 Double_t pt = d0tokpi->Pt();
398 Double_t rapidity = d0tokpi->YD0();
400 Double_t cosThetaStar = 9999.;
403 Double_t dca = d0tokpi->GetDCA();
406 Double_t d0xd0 = d0tokpi->Prodd0d0();
407 Double_t cosPointingAngle = d0tokpi->CosPointingAngle();
408 Double_t phi = d0tokpi->Phi();
409 Int_t pdgCode = mcVtxHF->GetPdgCode();
411 cosThetaStar = d0tokpi->CosThetaStarD0();
412 pTpi = d0tokpi->PtProng(0);
413 pTK = d0tokpi->PtProng(1);
414 d0pi = d0tokpi->Getd0Prong(0);
415 d0K = d0tokpi->Getd0Prong(1);
416 invMass=d0tokpi->InvMassD0();
419 cosThetaStar = d0tokpi->CosThetaStarD0bar();
420 pTpi = d0tokpi->PtProng(1);
421 pTK = d0tokpi->PtProng(0);
422 d0pi = d0tokpi->Getd0Prong(1);
423 d0K = d0tokpi->Getd0Prong(0);
424 invMass=d0tokpi->InvMassD0bar();
427 Double_t cT = d0tokpi->CtD0();
429 if (!fFillFromGenerated){
430 // ...either with reconstructed values....
431 containerInput[0] = pt;
432 containerInput[1] = rapidity;
433 containerInput[2] = cosThetaStar;
434 containerInput[3] = pTpi;
435 containerInput[4] = pTK;
436 containerInput[5] = cT*1.E4; // in micron
437 containerInput[6] = dca*1.E4; // in micron
438 containerInput[7] = d0pi*1.E4; // in micron
439 containerInput[8] = d0K*1.E4; // in micron
440 containerInput[9] = d0xd0*1.E8; // in micron^2
441 containerInput[10] = cosPointingAngle; // in micron
442 containerInput[11] = phi;
445 // ... or with generated values
446 Double_t vectorMC[7] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.};
447 if (GetGeneratedValuesFromMCParticle(mcVtxHF, mcArray, vectorMC)){
448 containerInput[0] = vectorMC[0];
449 containerInput[1] = vectorMC[1] ;
450 containerInput[2] = vectorMC[2] ;
451 containerInput[3] = vectorMC[3] ;
452 containerInput[4] = vectorMC[4] ;
453 containerInput[5] = vectorMC[5] ; // in micron
454 containerInput[6] = 0.; // dummy value, meaningless in MC, in micron
455 containerInput[7] = 0.; // dummy value, meaningless in MC, in micron
456 containerInput[8] = 0.; // dummy value, meaningless in MC, in micron
457 containerInput[9] = 100000.; // dummy value, meaningless in MC, in micron^2
458 containerInput[10] = 1.01; // dummy value, meaningless in MC
459 containerInput[11] = vectorMC[6];
462 AliDebug(3,"Problems in filling the container");
466 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]));
468 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
471 Bool_t acceptanceProng0 = (TMath::Abs(d0tokpi->EtaProng(0))< 0.9 && d0tokpi->PtProng(0) > 0.1);
472 Bool_t acceptanceProng1 = (TMath::Abs(d0tokpi->EtaProng(1))< 0.9 && d0tokpi->PtProng(1) > 0.1);
473 if (acceptanceProng0 && acceptanceProng1) {
474 AliDebug(2,"D0 reco daughters in acceptance");
475 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance) ;
480 Double_t fill[4]; //fill response matrix
482 // dimensions 0&1 : pt,eta (Rec)
487 // dimensions 2&3 : pt,eta (MC)
489 fill[2] = mcVtxHF->Pt();
490 fill[3] = mcVtxHF->Y();
492 fCorrelation->Fill(fill);
496 // cut on the min n. of clusters in ITS
497 Int_t ncls0=0,ncls1=0;
498 for(Int_t l=0;l<6;l++) {
499 if(TESTBIT(track0->GetITSClusterMap(),l)) ncls0++;
500 if(TESTBIT(track1->GetITSClusterMap(),l)) ncls1++;
502 AliDebug(2, Form("n clusters = %d", ncls0));
503 if (ncls0 >= fMinITSClusters && ncls1 >= fMinITSClusters) {
504 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
505 icountRecoITSClusters++;
506 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));
509 //added mass cut for D0 analysis
510 Double_t cuts[7] = {9999999., 1.1, 0., 9999999., 9999999., 0.,0.027};
511 //cuts of the D0 analysis (looser) see AliAnalysisTaskSED0Mass.cxx
515 cuts[0] = 400; //dca (um)
516 cuts[1] = 0.8; //cosTstar
517 cuts[2] = 0.5; //pt (GeV/c)
518 cuts[3] = 500; //d0 (um)
519 cuts[4] = -25000; //d0xd0 (um^2)
520 cuts[5] = 0.7; //cosTpointing
523 cuts[0] = 400; //dca (um)
524 cuts[1] = 0.8; //cosTstar
525 cuts[2] = 0.5; //pt (GeV/c)
526 cuts[3] = 500; //d0 (um)
527 cuts[4] = -20000; //d0xd0 (um^2)
528 cuts[5] = 0.5; //cosTpointing
531 /* //not same cuts for pt = 1 to 2 and 1 to 3 GeV in case must match them because of poor stat
532 else if (pt > 1 && pt <= 2){
552 else if (pt > 1 && pt <= 3){
571 else if (pt > 3 && pt <= 5){
610 Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
611 if (TMath::Abs(invMass-mD0PDG) < cuts[6]
613 && TMath::Abs(cosThetaStar) < cuts[1]
616 && TMath::Abs(d0pi*1E4) < cuts[3]
617 && TMath::Abs(d0K*1E4) < cuts[3]
618 && d0xd0*1E8 < cuts[4]
619 && cosPointingAngle > cuts[5]
622 AliDebug(2,"Particle passed PPR cuts (actually cuts for D0 analysis!)");
623 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPPR) ;
626 if(!fAcceptanceUnf){ // unfolding
628 Double_t fill[4]; //fill response matrix
630 // dimensions 0&1 : pt,eta (Rec)
635 // dimensions 2&3 : pt,eta (MC)
637 fill[2] = mcVtxHF->Pt();
638 fill[3] = mcVtxHF->Y();
640 fCorrelation->Fill(fill);
645 AliDebug(2,"Particle skipped due to PPR cuts");
646 if (dca*1E4 > cuts[0]){
647 AliDebug(2,"Problems with dca");
649 if ( TMath::Abs(cosThetaStar) > cuts[1]){
650 AliDebug(2,"Problems with cosThetaStar");
653 AliDebug(2,"Problems with pTpi");
656 AliDebug(2,"Problems with pTK");
658 if (TMath::Abs(d0pi*1E4) > cuts[3]){
659 AliDebug(2,"Problems with d0pi");
661 if (TMath::Abs(d0K*1E4) > cuts[3]){
662 AliDebug(2,"Problems with d0K");
664 if (d0xd0*1E8 > cuts[4]){
665 AliDebug(2,"Problems with d0xd0");
667 if (cosPointingAngle < cuts[5]){
668 AliDebug(2,"Problems with cosPointingAngle");
674 if(unsetvtx) d0tokpi->UnsetOwnPrimaryVtx();
675 } // end loop on D0->Kpi
677 AliDebug(2, Form("Found %i Reco particles that are D0!!",icountReco));
679 fCountReco+= icountReco;
680 fCountVertex+= icountVertex;
681 fCountRefit+= icountRefit;
682 fCountRecoAcc+= icountRecoAcc;
683 fCountRecoITSClusters+= icountRecoITSClusters;
684 fCountRecoPPR+= icountRecoPPR;
686 fHistEventsProcessed->Fill(0);
687 /* PostData(0) is taken care of by AliAnalysisTaskSE */
688 PostData(1,fHistEventsProcessed) ;
689 PostData(2,fCFManager->GetParticleContainer()) ;
690 PostData(3,fCorrelation) ;
694 //___________________________________________________________________________
695 void AliCFHeavyFlavourTaskMultiVarMultiStep::Terminate(Option_t*)
697 // The Terminate() function is the last function to be called during
698 // a query. It always runs on the client, it can be used to present
699 // the results graphically or save the results to file.
701 AliAnalysisTaskSE::Terminate();
703 AliInfo(Form("Found %i MC particles that are D0 in MC, in %d events",fCountMC,fEvents));
704 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, in %d events",fCountAcc,fEvents));
705 AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fEvents));
706 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));
707 AliInfo(Form("Found %i reco D0 that are decaying in K+pi, in %d events",fCountReco,fEvents));
708 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));
709 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));
710 AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and satisfy PPR cuts, in %d events",fCountRecoPPR,fEvents));
712 // draw some example plots....
714 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
716 printf("CONTAINER NOT FOUND\n");
719 // projecting the containers to obtain histograms
720 // first argument = variable, second argument = step
723 TH1D* h00 = cont->ShowProjection(0,0) ; // pt
724 TH1D* h10 = cont->ShowProjection(1,0) ; // rapidity
725 TH1D* h20 = cont->ShowProjection(2,0) ; // cosThetaStar
726 TH1D* h30 = cont->ShowProjection(3,0) ; // pTpi
727 TH1D* h40 = cont->ShowProjection(4,0) ; // pTK
728 TH1D* h50 = cont->ShowProjection(5,0) ; // cT
729 TH1D* h60 = cont->ShowProjection(6,0) ; // dca
730 TH1D* h70 = cont->ShowProjection(7,0) ; // d0pi
731 TH1D* h80 = cont->ShowProjection(8,0) ; // d0K
732 TH1D* h90 = cont->ShowProjection(9,0) ; // d0xd0
733 TH1D* h100 = cont->ShowProjection(10,0) ; // cosPointingAngle
734 TH1D* h110 = cont->ShowProjection(11,0) ; // phi
736 // MC-Acceptance level
737 TH1D* h01 = cont->ShowProjection(0,1) ; // pt
738 TH1D* h11 = cont->ShowProjection(1,1) ; // rapidity
739 TH1D* h21 = cont->ShowProjection(2,1) ; // cosThetaStar
740 TH1D* h31 = cont->ShowProjection(3,1) ; // pTpi
741 TH1D* h41 = cont->ShowProjection(4,1) ; // pTK
742 TH1D* h51 = cont->ShowProjection(5,1) ; // cT
743 TH1D* h61 = cont->ShowProjection(6,1) ; // dca
744 TH1D* h71 = cont->ShowProjection(7,1) ; // d0pi
745 TH1D* h81 = cont->ShowProjection(8,1) ; // d0K
746 TH1D* h91 = cont->ShowProjection(9,1) ; // d0xd0
747 TH1D* h101 = cont->ShowProjection(10,1) ; // cosPointingAngle
748 TH1D* h111 = cont->ShowProjection(11,1) ; // phi
751 TH1D* h04 = cont->ShowProjection(0,4) ; // pt
752 TH1D* h14 = cont->ShowProjection(1,4) ; // rapidity
753 TH1D* h24 = cont->ShowProjection(2,4) ; // cosThetaStar
754 TH1D* h34 = cont->ShowProjection(3,4) ; // pTpi
755 TH1D* h44 = cont->ShowProjection(4,4) ; // pTK
756 TH1D* h54 = cont->ShowProjection(5,4) ; // cT
757 TH1D* h64 = cont->ShowProjection(6,4) ; // dca
758 TH1D* h74 = cont->ShowProjection(7,4) ; // d0pi
759 TH1D* h84 = cont->ShowProjection(8,4) ; // d0K
760 TH1D* h94 = cont->ShowProjection(9,4) ; // d0xd0
761 TH1D* h104 = cont->ShowProjection(10,4) ; // cosPointingAngle
762 TH1D* h114 = cont->ShowProjection(11,4) ; // phi
764 h00->SetTitle("pT_D0 (GeV/c)");
765 h10->SetTitle("rapidity");
766 h20->SetTitle("cosThetaStar");
767 h30->SetTitle("pT_pi (GeV/c)");
768 h40->SetTitle("pT_K (Gev/c)");
769 h50->SetTitle("cT (#mum)");
770 h60->SetTitle("dca (#mum)");
771 h70->SetTitle("d0_pi (#mum)");
772 h80->SetTitle("d0_K (#mum)");
773 h90->SetTitle("d0xd0 (#mum^2)");
774 h100->SetTitle("cosPointingAngle");
775 h100->SetTitle("phi (rad)");
777 h00->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
778 h10->GetXaxis()->SetTitle("rapidity");
779 h20->GetXaxis()->SetTitle("cosThetaStar");
780 h30->GetXaxis()->SetTitle("pT_pi (GeV/c)");
781 h40->GetXaxis()->SetTitle("pT_K (Gev/c)");
782 h50->GetXaxis()->SetTitle("cT (#mum)");
783 h60->GetXaxis()->SetTitle("dca (#mum)");
784 h70->GetXaxis()->SetTitle("d0_pi (#mum)");
785 h80->GetXaxis()->SetTitle("d0_K (#mum)");
786 h90->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
787 h100->GetXaxis()->SetTitle("cosPointingAngle");
788 h110->GetXaxis()->SetTitle("phi (rad)");
790 h01->SetTitle("pT_D0 (GeV/c)");
791 h11->SetTitle("rapidity");
792 h21->SetTitle("cosThetaStar");
793 h31->SetTitle("pT_pi (GeV/c)");
794 h41->SetTitle("pT_K (Gev/c)");
795 h51->SetTitle("cT (#mum)");
796 h61->SetTitle("dca (#mum)");
797 h71->SetTitle("d0_pi (#mum)");
798 h81->SetTitle("d0_K (#mum)");
799 h91->SetTitle("d0xd0 (#mum^2)");
800 h101->SetTitle("cosPointingAngle");
801 h111->GetXaxis()->SetTitle("phi (rad)");
803 h01->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
804 h11->GetXaxis()->SetTitle("rapidity");
805 h21->GetXaxis()->SetTitle("cosThetaStar");
806 h31->GetXaxis()->SetTitle("pT_pi (GeV/c)");
807 h41->GetXaxis()->SetTitle("pT_K (Gev/c)");
808 h51->GetXaxis()->SetTitle("cT (#mum)");
809 h61->GetXaxis()->SetTitle("dca (#mum)");
810 h71->GetXaxis()->SetTitle("d0_pi (#mum)");
811 h81->GetXaxis()->SetTitle("d0_K (#mum)");
812 h91->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
813 h101->GetXaxis()->SetTitle("cosPointingAngle");
814 h111->GetXaxis()->SetTitle("phi (rad)");
816 h04->SetTitle("pT_D0 (GeV/c)");
817 h14->SetTitle("rapidity");
818 h24->SetTitle("cosThetaStar");
819 h34->SetTitle("pT_pi (GeV/c)");
820 h44->SetTitle("pT_K (Gev/c)");
821 h54->SetTitle("cT (#mum)");
822 h64->SetTitle("dca (#mum)");
823 h74->SetTitle("d0_pi (#mum)");
824 h84->SetTitle("d0_K (#mum)");
825 h94->SetTitle("d0xd0 (#mum^2)");
826 h104->SetTitle("cosPointingAngle");
827 h114->GetXaxis()->SetTitle("phi (rad)");
829 h04->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
830 h14->GetXaxis()->SetTitle("rapidity");
831 h24->GetXaxis()->SetTitle("cosThetaStar");
832 h34->GetXaxis()->SetTitle("pT_pi (GeV/c)");
833 h44->GetXaxis()->SetTitle("pT_K (Gev/c)");
834 h54->GetXaxis()->SetTitle("cT (#mum)");
835 h64->GetXaxis()->SetTitle("dca (#mum)");
836 h74->GetXaxis()->SetTitle("d0_pi (#mum)");
837 h84->GetXaxis()->SetTitle("d0_K (#mum)");
838 h94->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
839 h104->GetXaxis()->SetTitle("cosPointingAngle");
840 h114->GetXaxis()->SetTitle("phi (rad)");
842 Double_t max0 = h00->GetMaximum();
843 Double_t max1 = h10->GetMaximum();
844 Double_t max2 = h20->GetMaximum();
845 Double_t max3 = h30->GetMaximum();
846 Double_t max4 = h40->GetMaximum();
847 Double_t max5 = h50->GetMaximum();
848 Double_t max6 = h60->GetMaximum();
849 Double_t max7 = h70->GetMaximum();
850 Double_t max8 = h80->GetMaximum();
851 Double_t max9 = h90->GetMaximum();
852 Double_t max10 = h100->GetMaximum();
853 Double_t max11 = h110->GetMaximum();
855 h00->GetYaxis()->SetRangeUser(0,max0*1.2);
856 h10->GetYaxis()->SetRangeUser(0,max1*1.2);
857 h20->GetYaxis()->SetRangeUser(0,max2*1.2);
858 h30->GetYaxis()->SetRangeUser(0,max3*1.2);
859 h40->GetYaxis()->SetRangeUser(0,max4*1.2);
860 h50->GetYaxis()->SetRangeUser(0,max5*1.2);
861 h60->GetYaxis()->SetRangeUser(0,max6*1.2);
862 h70->GetYaxis()->SetRangeUser(0,max7*1.2);
863 h80->GetYaxis()->SetRangeUser(0,max8*1.2);
864 h90->GetYaxis()->SetRangeUser(0,max9*1.2);
865 h100->GetYaxis()->SetRangeUser(0,max10*1.2);
866 h110->GetYaxis()->SetRangeUser(0,max11*1.2);
868 h01->GetYaxis()->SetRangeUser(0,max0*1.2);
869 h11->GetYaxis()->SetRangeUser(0,max1*1.2);
870 h21->GetYaxis()->SetRangeUser(0,max2*1.2);
871 h31->GetYaxis()->SetRangeUser(0,max3*1.2);
872 h41->GetYaxis()->SetRangeUser(0,max4*1.2);
873 h51->GetYaxis()->SetRangeUser(0,max5*1.2);
874 h61->GetYaxis()->SetRangeUser(0,max6*1.2);
875 h71->GetYaxis()->SetRangeUser(0,max7*1.2);
876 h81->GetYaxis()->SetRangeUser(0,max8*1.2);
877 h91->GetYaxis()->SetRangeUser(0,max9*1.2);
878 h101->GetYaxis()->SetRangeUser(0,max10*1.2);
879 h111->GetYaxis()->SetRangeUser(0,max11*1.2);
881 h00->SetMarkerStyle(20);
882 h10->SetMarkerStyle(24);
883 h20->SetMarkerStyle(21);
884 h30->SetMarkerStyle(25);
885 h40->SetMarkerStyle(27);
886 h50->SetMarkerStyle(28);
887 h60->SetMarkerStyle(20);
888 h70->SetMarkerStyle(24);
889 h80->SetMarkerStyle(21);
890 h90->SetMarkerStyle(25);
891 h100->SetMarkerStyle(27);
892 h110->SetMarkerStyle(28);
894 h00->SetMarkerColor(2);
895 h10->SetMarkerColor(2);
896 h20->SetMarkerColor(2);
897 h30->SetMarkerColor(2);
898 h40->SetMarkerColor(2);
899 h50->SetMarkerColor(2);
900 h60->SetMarkerColor(2);
901 h70->SetMarkerColor(2);
902 h80->SetMarkerColor(2);
903 h90->SetMarkerColor(2);
904 h100->SetMarkerColor(2);
905 h110->SetMarkerColor(2);
907 h01->SetMarkerStyle(20) ;
908 h11->SetMarkerStyle(24) ;
909 h21->SetMarkerStyle(21) ;
910 h31->SetMarkerStyle(25) ;
911 h41->SetMarkerStyle(27) ;
912 h51->SetMarkerStyle(28) ;
913 h61->SetMarkerStyle(20);
914 h71->SetMarkerStyle(24);
915 h81->SetMarkerStyle(21);
916 h91->SetMarkerStyle(25);
917 h101->SetMarkerStyle(27);
918 h111->SetMarkerStyle(28);
920 h01->SetMarkerColor(8);
921 h11->SetMarkerColor(8);
922 h21->SetMarkerColor(8);
923 h31->SetMarkerColor(8);
924 h41->SetMarkerColor(8);
925 h51->SetMarkerColor(8);
926 h61->SetMarkerColor(8);
927 h71->SetMarkerColor(8);
928 h81->SetMarkerColor(8);
929 h91->SetMarkerColor(8);
930 h101->SetMarkerColor(8);
931 h111->SetMarkerColor(8);
933 h04->SetMarkerStyle(20) ;
934 h14->SetMarkerStyle(24) ;
935 h24->SetMarkerStyle(21) ;
936 h34->SetMarkerStyle(25) ;
937 h44->SetMarkerStyle(27) ;
938 h54->SetMarkerStyle(28) ;
939 h64->SetMarkerStyle(20);
940 h74->SetMarkerStyle(24);
941 h84->SetMarkerStyle(21);
942 h94->SetMarkerStyle(25);
943 h104->SetMarkerStyle(27);
944 h114->SetMarkerStyle(28);
946 h04->SetMarkerColor(4);
947 h14->SetMarkerColor(4);
948 h24->SetMarkerColor(4);
949 h34->SetMarkerColor(4);
950 h44->SetMarkerColor(4);
951 h54->SetMarkerColor(4);
952 h64->SetMarkerColor(4);
953 h74->SetMarkerColor(4);
954 h84->SetMarkerColor(4);
955 h94->SetMarkerColor(4);
956 h104->SetMarkerColor(4);
957 h114->SetMarkerColor(4);
959 gStyle->SetCanvasColor(0);
960 gStyle->SetFrameFillColor(0);
961 gStyle->SetTitleFillColor(0);
962 gStyle->SetStatColor(0);
964 // drawing in 2 separate canvas for a matter of clearity
965 TCanvas * c1 =new TCanvas("c1","pT, rapidiy, cosThetaStar",1100,1600);
997 TCanvas * c2 =new TCanvas("c2","pTpi, pTK, cT",1100,1600);
1028 TCanvas * c3 =new TCanvas("c3","dca, d0pi, d0K",1100,1600);
1059 TCanvas * c4 =new TCanvas("c4","d0xd0, cosPointingAngle, phi",1100,1600);
1090 THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
1092 TH2D* corr1 =hcorr->Projection(0,2);
1093 TH2D* corr2 = hcorr->Projection(1,3);
1095 TCanvas * c7 =new TCanvas("c7","",800,400);
1098 corr1->Draw("text");
1100 corr2->Draw("text");
1103 TFile* file_projection = new TFile("file_projection.root","RECREATE");
1107 h00->Write("pT_D0_step0");
1108 h10->Write("rapidity_step0");
1109 h20->Write("cosThetaStar_step0");
1110 h30->Write("pT_pi_step0");
1111 h40->Write("pT_K_step0");
1112 h50->Write("cT_step0");
1113 h60->Write("dca_step0");
1114 h70->Write("d0_pi_step0");
1115 h80->Write("d0_K_step0");
1116 h90->Write("d0xd0_step0");
1117 h100->Write("cosPointingAngle_step0");
1118 h110->Write("phi_step0");
1120 h01->Write("pT_D0_step1");
1121 h11->Write("rapidity_step1");
1122 h21->Write("cosThetaStar_step1");
1123 h31->Write("pT_pi_step1");
1124 h41->Write("pT_K_step1");
1125 h51->Write("cT_step1");
1126 h61->Write("dca_step1");
1127 h71->Write("d0_pi_step1");
1128 h81->Write("d0_K_step1");
1129 h91->Write("d0xd0_step1");
1130 h101->Write("cosPointingAngle_step1");
1131 h111->Write("phi_step1");
1133 h04->Write("pT_D0_step2");
1134 h14->Write("rapidity_step2");
1135 h24->Write("cosThetaStar_step2");
1136 h34->Write("pT_pi_step2");
1137 h44->Write("pT_K_step2");
1138 h54->Write("cT_step2");
1139 h64->Write("dca_step2");
1140 h74->Write("d0_pi_step2");
1141 h80->Write("d0_K_step2");
1142 h94->Write("d0xd0_step2");
1143 h104->Write("cosPointingAngle_step2");
1144 h114->Write("phi_step2");
1146 file_projection->Close();
1151 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
1152 c2->SaveAs("Plots/pTpi_pTK_cT.eps");
1153 c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
1154 c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
1156 c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
1157 c2->SaveAs("Plots/pTpi_pTK_cT.gif");
1158 c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
1159 c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
1163 //___________________________________________________________________________
1164 void AliCFHeavyFlavourTaskMultiVarMultiStep::UserCreateOutputObjects() {
1165 //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1166 //TO BE SET BEFORE THE EXECUTION OF THE TASK
1168 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1172 fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
1175 //___________________________________________________________________________
1176 Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CosThetaStar(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
1179 // to calculate cos(ThetaStar) for generated particle
1180 // using the K, since mcPartDaughter0 and mcPartDaughter1 always correspond to K and pi respectively
1181 // (see where the function is called)
1184 Int_t pdgvtx = mcPart->GetPdgCode();
1185 /* if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1186 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1187 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1188 AliInfo(Form("D0, with pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1189 AliDebug(2,"This is a D0");
1190 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1191 mcPartDaughter0 = mcPartDaughter1;
1192 mcPartDaughter1 = mcPartdummy;
1198 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1199 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1200 if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1204 AliDebug(2,"D0bar");
1206 if (pdgprong0 == 211){
1207 AliDebug(2,Form("pdgprong0 = %d, pdgprong1 = %d, switching...",pdgprong0,pdgprong1));
1208 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1209 mcPartDaughter0 = mcPartDaughter1;
1210 mcPartDaughter1 = mcPartdummy;
1211 pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1212 pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1215 AliDebug(2,Form("After checking, pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1216 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdgvtx))->Mass();
1218 massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
1219 massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
1221 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);
1223 Double_t px = mcPartDaughter0->Px()+mcPartDaughter1->Px();
1224 Double_t py = mcPartDaughter0->Py()+mcPartDaughter1->Py();
1225 Double_t pz = mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
1226 Double_t p = TMath::Sqrt(px*px+py*py+pz*pz);
1227 Double_t e = TMath::Sqrt(massvtx*massvtx+p*p);
1228 Double_t beta = p/e;
1229 Double_t gamma = e/massvtx;
1230 TVector3 mom(mcPartDaughter0->Px(),mcPartDaughter0->Py(),mcPartDaughter0->Pz());
1231 TVector3 momTot(mcPartDaughter0->Px()+mcPartDaughter1->Px(),mcPartDaughter0->Py()+mcPartDaughter1->Py(),mcPartDaughter0->Pz()+mcPartDaughter1->Pz());
1233 Double_t qlprong = mom.Dot(momTot)/momTot.Mag(); // analog to AliAODRecoDecay::QlProng(0)
1235 AliDebug(2,Form("pStar = %f, beta = %f, gamma = %f, qlprong = %f, massp[0] = %f", pStar, beta, gamma, qlprong, massp[0]));
1236 Double_t cts = (qlprong/gamma-beta*TMath::Sqrt(pStar*pStar+massp[0]*massp[0]))/pStar;
1237 AliDebug(2,Form("cts = %f", cts));
1240 //___________________________________________________________________________
1241 Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CT(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
1244 // to calculate cT for generated particle
1247 Double_t xmcPart[3] = {0,0,0};
1248 Double_t xdaughter0[3] = {0,0,0};
1249 Double_t xdaughter1[3] = {0,0,0};
1250 mcPart->XvYvZv(xmcPart); // cm
1251 mcPartDaughter0->XvYvZv(xdaughter0); // cm
1252 mcPartDaughter1->XvYvZv(xdaughter1); //cm
1253 Double_t prodVtxD0 = TMath::Sqrt(xmcPart[0]*xmcPart[0]+
1254 xmcPart[1]*xmcPart[1]+
1255 xmcPart[2]*xmcPart[2]);
1256 Double_t prodVtxDaughter0 = TMath::Sqrt(xdaughter0[0]*xdaughter0[0]+
1257 xdaughter0[1]*xdaughter0[1]+
1258 xdaughter0[2]*xdaughter0[2]);
1259 Double_t prodVtxDaughter1 = TMath::Sqrt(xdaughter1[0]*xdaughter1[0]+
1260 xdaughter1[1]*xdaughter1[1]+
1261 xdaughter1[2]*xdaughter1[2]);
1263 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));
1264 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));
1265 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));
1267 Double_t cT0 = TMath::Sqrt((xdaughter0[0]-xmcPart[0])*(xdaughter0[0]-xmcPart[0])+
1268 (xdaughter0[1]-xmcPart[1])*(xdaughter0[1]-xmcPart[1])+
1269 (xdaughter0[2]-xmcPart[2])*(xdaughter0[2]-xmcPart[2]));
1271 Double_t cT1 = TMath::Sqrt((xdaughter1[0]-xmcPart[0])*(xdaughter1[0]-xmcPart[0])+
1272 (xdaughter1[1]-xmcPart[1])*(xdaughter1[1]-xmcPart[1])+
1273 (xdaughter1[2]-xmcPart[2])*(xdaughter1[2]-xmcPart[2]));
1275 if ((cT0 - cT1)>1E-5) {
1276 AliWarning(Form("cT from daughter 0 (%f) different from cT from daughter 1 (%f)! Using cT from daughter 0, but PLEASE, CHECK....",cT0,cT1));
1279 // calculating cT from cT0
1281 Double_t mass = TDatabasePDG::Instance()->GetParticle(mcPart->GetPdgCode())->Mass(); // mcPart->GetPdgCode() should return +/- 421 for the D0/D0bar
1282 Double_t p = mcPart-> P();
1283 Double_t cT = cT0*mass/p;
1284 AliDebug(2, Form("cT from daughter 0 = %f (micron)", cT0*1E4));
1285 AliDebug(2, Form("cT from daughter 1 = %f (micron)", cT1*1E4));
1286 AliDebug(2, Form("cT (with mass = %f and p = %f) = %f (micron)", mass, p, cT*1E4));
1289 //________________________________________________________________________________
1290 Bool_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC) const {
1293 // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
1296 Bool_t isOk = kFALSE;
1298 // check whether the D0 decays in pi+K
1299 // to be added!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1300 // could use a cut in the CF, but not clear how to define a TDecayChannel
1301 // implemented for the time being as a cut on the number of daughters - see later when
1302 // getting the daughters
1304 // getting the daughters
1305 Int_t daughter0 = mcPart->GetDaughter(0);
1306 Int_t daughter1 = mcPart->GetDaughter(1);
1307 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
1308 if (daughter0 == 0 || daughter1 == 0) {
1309 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
1312 if (TMath::Abs(daughter1 - daughter0) != 1) {
1313 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
1316 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
1317 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
1318 if (!mcPartDaughter0 || !mcPartDaughter1) {
1319 AliWarning("At least one Daughter Particle not found in tree, skipping");
1322 if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==321 &&
1323 TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
1324 !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
1325 TMath::Abs(mcPartDaughter1->GetPdgCode())==321)) {
1326 AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
1330 Double_t vtx1[3] = {0,0,0}; // primary vertex
1331 Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
1332 Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
1333 mcPart->XvYvZv(vtx1); // cm
1334 // getting vertex from daughters
1335 mcPartDaughter0->XvYvZv(vtx2daughter0); // cm
1336 mcPartDaughter1->XvYvZv(vtx2daughter1); //cm
1337 if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
1338 AliError("Daughters have different secondary vertex, skipping the track");
1343 // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
1344 AliAODMCParticle* positiveDaugh = mcPartDaughter0;
1345 AliAODMCParticle* negativeDaugh = mcPartDaughter1;
1346 if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
1347 // inverting in case the positive daughter is the second one
1348 positiveDaugh = mcPartDaughter1;
1349 negativeDaugh = mcPartDaughter0;
1351 // getting the momentum from the daughters
1352 Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
1353 Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
1354 Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
1356 Double_t d0[2] = {0.,0.};
1358 AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
1360 Double_t cosThetaStar = 0.;
1361 Double_t cosThetaStarD0 = 0.;
1362 Double_t cosThetaStarD0bar = 0.;
1363 cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
1364 cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
1365 if (mcPart->GetPdgCode() == 421){ // D0
1366 AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1367 cosThetaStar = cosThetaStarD0;
1369 else if (mcPart->GetPdgCode() == -421){ // D0bar{
1370 AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1371 cosThetaStar = cosThetaStarD0bar;
1374 AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
1377 if (cosThetaStar < -1 || cosThetaStar > 1) {
1378 AliWarning("Invalid value for cosine Theta star, returning");
1382 // calculate cos(Theta*) according to the method implemented herein
1384 Double_t cts = 9999.;
1385 cts = CosThetaStar(mcPart, mcPartDaughter0, mcPartDaughter1);
1386 if (cts < -1 || cts > 1) {
1387 AliWarning("Invalid value for cosine Theta star from AliCFHeavyFlavourTaskMultiVarMultiStep method");
1389 if (TMath::Abs(cts - cosThetaStar)>0.001) {
1390 AliError(Form("cosThetaStar automatically calculated different from that manually calculated!!! cosThetaStar = %f, cosThetaStar = %f", cosThetaStar,cts));
1393 Double_t cT = decay->Ct(421);
1395 // calculate cT from the method implemented herein
1397 cT1 = CT(mcPart, mcPartDaughter0, mcPartDaughter1);
1399 if (TMath::Abs(cT1 - cT)>0.001) {
1400 AliError(Form("cT automatically calculated different from that manually calculated!!! cT = %f, cT1 = %f",cT,cT1));
1403 // get the pT of the daughters
1408 if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
1409 pTpi = mcPartDaughter0->Pt();
1410 pTK = mcPartDaughter1->Pt();
1413 pTpi = mcPartDaughter1->Pt();
1414 pTK = mcPartDaughter0->Pt();
1417 vectorMC[0] = mcPart->Pt();
1418 vectorMC[1] = mcPart->Y() ;
1419 vectorMC[2] = cosThetaStar ;
1420 vectorMC[3] = pTpi ;
1422 vectorMC[5] = cT*1.E4 ; // in micron
1423 vectorMC[6] = mcPart->Phi() ;
1427 //_________________________________________________________________________________________________
1428 Int_t AliCFHeavyFlavourTaskMultiVarMultiStep::CheckOrigin(AliAODMCParticle* mcPart, TClonesArray* mcArray)const{
1431 // checking whether the very mother of the D0 is a charm or a bottom quark
1434 Int_t pdgGranma = 0;
1436 mother = mcPart->GetMother();
1440 AliDebug(2,Form("mother at step %d = %d", istep, mother));
1441 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
1442 pdgGranma = mcGranma->GetPdgCode();
1443 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1444 Int_t abspdgGranma = TMath::Abs(pdgGranma);
1445 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1448 mother = mcGranma->GetMother();