1 /**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 //-----------------------------------------------------------------------
19 // Class for HF corrections as a function of many variables and step
20 // Author : C. Zampolli, CERN
21 // D. Caffarri, Univ & INFN Padova caffarri@pd.infn.it
22 // Base class for HF Unfolding - agrelli@uu.nl
23 //-----------------------------------------------------------------------
25 #include "TParticle.h"
26 #include "TClonesArray.h"
27 #include "AliAODMCParticle.h"
28 #include "AliAODRecoDecayHF.h"
29 #include "AliAODRecoDecayHF2Prong.h"
30 #include "AliAODRecoDecayHF3Prong.h"
31 #include "AliAODRecoDecayHF4Prong.h"
32 #include "AliAODRecoCascadeHF.h"
33 #include "AliAODMCHeader.h"
34 #include "AliAODEvent.h"
36 #include "AliESDtrackCuts.h"
37 #include "AliESDtrack.h"
39 #include "AliCFVertexingHF.h"
41 //___________________________________________________________
42 AliCFVertexingHF::AliCFVertexingHF() :
45 fmcPartCandidate(0x0),
50 fFillFromGenerated(0),
53 fKeepDfromBOnly(kFALSE),
61 fFake(1.), // setting to MC value
62 fRejectIfNoQuark(kFALSE)
74 //_____________________________________________________
75 AliCFVertexingHF::AliCFVertexingHF(TClonesArray *mcArray, UShort_t originDselection) :
78 fmcPartCandidate(0x0),
83 fFillFromGenerated(0),
86 fKeepDfromBOnly(kFALSE),
94 fFake(1.), // setting to MC value
95 fRejectIfNoQuark(kFALSE)
98 // constructor with mcArray
101 SetDselection(originDselection);
105 //_______________________________________________________
106 AliCFVertexingHF::~AliCFVertexingHF()
112 if (fmcArray) fmcArray = 0x0;
113 if (fRecoCandidate) fRecoCandidate = 0x0;
114 if (fmcPartCandidate) fmcPartCandidate = 0x0;
116 delete [] fLabelArray;
124 delete [] fEtaAccCut;
129 //_____________________________________________________
130 AliCFVertexingHF& AliCFVertexingHF::operator=(const AliCFVertexingHF& c)
133 // assigment operator
137 TObject::operator=(c);
138 fmcArray = c.fmcArray;
139 fRecoCandidate = c.fRecoCandidate;
140 fmcPartCandidate = c.fmcPartCandidate;
141 fNDaughters = c.fNDaughters;
143 fzPrimVertex = c.fzPrimVertex;
144 fzMCVertex = c.fzMCVertex;
145 fFillFromGenerated = c.fFillFromGenerated;
146 fOriginDselection = c.fOriginDselection;
147 fKeepDfromB = c.fKeepDfromB;
148 fKeepDfromBOnly = c.fKeepDfromBOnly;
149 fmcLabel = c.fmcLabel;
151 fCentValue=c.fCentValue;
152 fFakeSelection=c.fFakeSelection;
154 fRejectIfNoQuark=c.fRejectIfNoQuark;
156 fLabelArray = new Int_t[fProngs];
157 fPtAccCut = new Float_t[fProngs];
158 fEtaAccCut = new Float_t[fProngs];
159 for(Int_t iP=0; iP<fProngs; iP++){
160 fLabelArray[iP]=c.fLabelArray[iP];
161 fPtAccCut[iP]=c.fPtAccCut[iP];
162 fEtaAccCut[iP]=c.fEtaAccCut[iP];
170 //____________________________________________________
171 AliCFVertexingHF::AliCFVertexingHF(const AliCFVertexingHF &c) :
173 fmcArray(c.fmcArray),
174 fRecoCandidate(c.fRecoCandidate),
175 fmcPartCandidate(c.fmcPartCandidate),
176 fNDaughters(c.fNDaughters),
178 fzPrimVertex(c.fzPrimVertex),
179 fzMCVertex(c.fzMCVertex),
180 fFillFromGenerated(c.fFillFromGenerated),
181 fOriginDselection (c.fOriginDselection),
182 fKeepDfromB (c.fKeepDfromB),
183 fKeepDfromBOnly (c.fKeepDfromBOnly),
184 fmcLabel(c.fmcLabel),
187 fCentValue(c.fCentValue),
190 fFakeSelection(c.fFakeSelection),
192 fRejectIfNoQuark(c.fRejectIfNoQuark)
198 fLabelArray = new Int_t[fProngs];
199 fPtAccCut = new Float_t[fProngs];
200 fEtaAccCut = new Float_t[fProngs];
201 if (c.fLabelArray) memcpy(fLabelArray,c.fLabelArray,fProngs*sizeof(Int_t));
202 if (c.fPtAccCut) memcpy(fPtAccCut,c.fPtAccCut,fProngs*sizeof(Int_t));
203 if (c.fEtaAccCut) memcpy(fEtaAccCut,c.fEtaAccCut,fProngs*sizeof(Int_t));
207 //___________________________________________________________
208 void AliCFVertexingHF::SetDselection(UShort_t originDselection)
210 // setting the way the D0 will be selected
211 // 0 --> only from c quarks
212 // 1 --> only from b quarks
213 // 2 --> from both c quarks and b quarks
215 fOriginDselection = originDselection;
217 if (fOriginDselection == 0) {
218 fKeepDfromB = kFALSE;
219 fKeepDfromBOnly = kFALSE;
222 if (fOriginDselection == 1) {
224 fKeepDfromBOnly = kTRUE;
227 if (fOriginDselection == 2) {
229 fKeepDfromBOnly = kFALSE;
235 //______________________________________________________
236 void AliCFVertexingHF::SetMCCandidateParam(Int_t label)
239 // setting the parameters (candidate and n. daughters)
242 fmcPartCandidate = dynamic_cast <AliAODMCParticle*> (fmcArray->At(label));
243 if (fmcPartCandidate){
244 fNDaughters = fmcPartCandidate->GetNDaughters();
247 AliError(Form("Dynamic cast failed, fNdaughters will remain set to %d",fNDaughters));
252 //____________________________________________________________
253 Int_t AliCFVertexingHF::MCcquarkCounting(AliAODMCParticle* mcPart) const
256 // counting the c-quarks
261 if (mcPart->GetPdgCode() == 4) cquarks++;
262 if (mcPart->GetPdgCode() == -4) cquarks++;
265 AliWarning("Particle not found in tree, skipping\n");
272 //________________________________________________________
273 Bool_t AliCFVertexingHF::CheckMCPartFamily(AliAODMCParticle */*mcPart*/, TClonesArray */*mcArray*/) const
276 //checking the family
279 Int_t pdgGranma = CheckOrigin();
281 if (pdgGranma == -99999){
282 AliDebug(2,"This particle does not have a quark in his genealogy\n");
285 if (pdgGranma == -9999){
286 AliDebug(2,"This particle come from a B decay channel but according to the settings of the task, we keep only the prompt charm particles\n");
290 if (pdgGranma == -999){
291 AliDebug(2,"This particle come from a prompt charm particles but according to the settings of the task, we want only the ones coming from B\n");
295 if (!CheckMCDaughters()) {
296 AliDebug(3, "CheckMCDaughters false");
299 if (!CheckMCChannelDecay()) {
300 AliDebug(3,"CheckMCChannelDecay false");
306 //_________________________________________________________________________________________________
307 Int_t AliCFVertexingHF::CheckOrigin() const
310 // checking whether the mother of the particles come from a charm or a bottom quark
315 mother = fmcPartCandidate->GetMother();
317 Int_t abspdgGranma =0;
318 Bool_t isFromB=kFALSE;
319 Bool_t isQuarkFound=kFALSE;
322 AliDebug(2,Form("mother at step %d = %d", istep, mother));
323 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(fmcArray->At(mother));
325 pdgGranma = mcGranma->GetPdgCode();
326 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
327 abspdgGranma = TMath::Abs(pdgGranma);
328 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
331 if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
332 mother = mcGranma->GetMother();
334 AliError("Failed casting the mother particle!");
339 if(fRejectIfNoQuark && !isQuarkFound) return -99999;
341 if (!fKeepDfromB) return -9999; //skip particle if come from a B meson.
344 if (fKeepDfromBOnly) return -999;
349 //___________________________________________
350 Bool_t AliCFVertexingHF::CheckMCDaughters()const
353 // checking the daughters
356 AliAODMCParticle *mcPartDaughter;
357 Bool_t checkDaughters = kFALSE;
359 Int_t label0 = fmcPartCandidate->GetDaughter(0);
360 Int_t label1 = fmcPartCandidate->GetDaughter(1);
361 AliDebug(3,Form("label0 = %d, label1 = %d",label0,label1));
362 if (label1==0 || label0 == 0){
363 AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
364 return checkDaughters;
367 if (fLabelArray == 0x0) {
368 return checkDaughters;
371 for (Int_t iProng = 0; iProng<fProngs; iProng++){
372 mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fLabelArray[iProng]));
373 if (!mcPartDaughter) {
374 AliWarning("At least one Daughter Particle not found in tree, skipping");
375 return checkDaughters;
379 checkDaughters = kTRUE;
380 return checkDaughters;
383 //______________________________________________________
384 Bool_t AliCFVertexingHF::FillMCContainer(Double_t *containerInputMC)
387 // fill the container for Generator level selection
390 Bool_t mcContainerFilled = kFALSE;
392 Double_t* vectorMC = new Double_t[fNVar];
393 for (Int_t iVar = 0; iVar<fNVar; iVar++) vectorMC[iVar]= 9999.;
395 if (GetGeneratedValuesFromMCParticle(&vectorMC[0])){
396 for (Int_t iVar = 0; iVar<fNVar; iVar++){
397 containerInputMC[iVar] = vectorMC[iVar];
399 mcContainerFilled = kTRUE;
403 return mcContainerFilled;
406 //______________________________________________________
407 Bool_t AliCFVertexingHF::FillRecoContainer(Double_t *containerInput)
410 // fill the container for Reconstrucred level selection
413 Bool_t recoContainerFilled = kFALSE;
414 Double_t* vectorValues = new Double_t[fNVar];
415 Double_t* vectorReco = new Double_t[fNVar];
416 for (Int_t iVar = 0; iVar<fNVar; iVar++) {
418 vectorValues[iVar]= 9999.;
419 vectorReco[iVar]=9999.;
422 if(fFillFromGenerated){
423 //filled with MC values
424 if (GetGeneratedValuesFromMCParticle(&vectorValues[0])){
425 for (Int_t iVar = 0; iVar<fNVar; iVar++){
426 containerInput[iVar] = vectorValues[iVar];
428 recoContainerFilled = kTRUE;
432 //filled with Reco values
434 if (GetRecoValuesFromCandidate(&vectorReco[0])){
435 for (Int_t iVar = 0; iVar<fNVar; iVar++){
436 containerInput[iVar] = vectorReco[iVar];
438 recoContainerFilled = kTRUE;
442 delete [] vectorValues;
443 delete [] vectorReco;
446 return recoContainerFilled;
449 //_____________________________________________________
450 Bool_t AliCFVertexingHF::MCAcceptanceStep() const
453 // checking the MC acceptance step
456 Bool_t bMCAccStep = kFALSE;
458 AliAODMCParticle *mcPartDaughter;
459 Int_t label0 = fmcPartCandidate->GetDaughter(0);
460 Int_t label1 = fmcPartCandidate->GetDaughter(1);
461 if (label1==0 || label0 == 0){
462 AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
466 if (fLabelArray == 0x0) {
470 for (Int_t iProng = 0; iProng<fProngs; iProng++){
471 mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(fLabelArray[iProng]));
472 if (!mcPartDaughter) {
473 AliWarning("At least one Daughter Particle not found in tree, skipping");
476 Double_t eta = mcPartDaughter->Eta();
477 Double_t pt = mcPartDaughter->Pt();
479 //set values of eta and pt in the constructor.
480 // if (TMath::Abs(eta) > 0.9 || pt < 0.1){
481 if (TMath::Abs(eta) > fEtaAccCut[iProng] || pt < fPtAccCut[iProng]){
482 AliDebug(3,Form("At least one daughter has eta or pt outside the required range (|eta| = %f, pt = %f, should be |eta| < %f, pt > %f \n", TMath::Abs(eta), pt, fEtaAccCut[iProng], fPtAccCut[iProng]));
490 //_____________________________________________________
491 Bool_t AliCFVertexingHF::MCRefitStep(AliAODEvent *aodEvent, AliESDtrackCuts **trackCuts) const
494 // check on the kTPCrefit and kITSrefit conditions of the daughters
496 Bool_t bRefitStep = kFALSE;
498 Int_t label0 = fmcPartCandidate->GetDaughter(0);
499 Int_t label1 = fmcPartCandidate->GetDaughter(1);
501 if (label1==0 || label0 == 0){
502 AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
506 if (fLabelArray == 0x0) {
510 Int_t foundDaughters = 0;
511 Int_t* temp = new Int_t[fProngs];
512 for (Int_t ilabel = 0; ilabel<fProngs; ilabel++){
513 temp[ilabel] = fLabelArray[ilabel];
516 // if (trackCuts->GetRequireTPCRefit() || trackCuts->GetRequireITSRefit()){
518 for(Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
519 AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
520 if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
521 Bool_t foundTrack = kFALSE;
523 for (Int_t ilabel = 0; ilabel<fProngs; ilabel++){
524 AliDebug(3,Form("fLabelArray[%d] = %d, track->GetLabel() = %d",ilabel,fLabelArray[ilabel],TMath::Abs(track->GetLabel())));
525 if ((track->GetLabel()<0)&&(fFakeSelection==1)) continue;
526 if ((track->GetLabel()>0)&&(fFakeSelection==2)) continue;
528 if (TMath::Abs(track->GetLabel()) == temp[ilabel]) {
537 AliDebug(4,Form("daughter %d \n",foundDaughters));
538 if(trackCuts[prongindex]->GetRequireTPCRefit()){
539 if(track->GetStatus()&AliESDtrack::kTPCrefit) {
543 AliDebug(3, "Refit cut not passed , missing TPC refit\n");
550 if (trackCuts[prongindex]->GetRequireITSRefit()) {
551 if(track->GetStatus()&AliESDtrack::kITSrefit){
555 AliDebug(3, "Refit cut not passed , missing ITS refit\n");
562 if (foundDaughters == fProngs){
569 if (foundDaughters== fProngs) return bRefitStep;
573 //____________________________________________________________________________
575 Bool_t AliCFVertexingHF::RecoStep()
578 //check also vertex and ITS Refit and TPC Refit
581 Bool_t bRecoStep = kFALSE;
582 Int_t mcLabel = GetMCLabel();
585 AliDebug(2,"No MC particle found");
589 fmcPartCandidate = (AliAODMCParticle*)fmcArray->At(mcLabel);
590 if (!fmcPartCandidate){
591 AliWarning("Could not find associated MC in AOD MC tree");
596 Int_t pdgGranma = CheckOrigin();
598 if (pdgGranma == -99999){
599 AliDebug(2,"This particle does not have a quark in his genealogy\n");
602 if (pdgGranma == -9999){
603 AliDebug(2,"This particle come from a B decay channel but according to the settings of the task, we keep only prompt charm particles\n");
607 if (pdgGranma == -999){
608 AliDebug(2,"This particle come from a prompt charm particle but according to the settings of the task, we want only the ones coming from B\n");
615 //____________________________________________
616 Double_t AliCFVertexingHF::GetEtaProng(Int_t iProng) const
619 // getting eta of the prong
623 Double_t etaProng = fRecoCandidate->EtaProng(iProng);
628 //______________________________________________________
629 Double_t AliCFVertexingHF::GetPtProng(Int_t iProng) const
632 // getting pt of the prong
636 Double_t ptProng = fRecoCandidate->PtProng(iProng);
643 //____________________________________________________________________
645 Bool_t AliCFVertexingHF::RecoAcceptStep(AliESDtrackCuts **trackCuts) const
648 // reco Acceptance step
651 Bool_t bRecoAccStep = kFALSE;
653 Float_t etaCutMin, ptCutMin, etaCutMax, ptCutMax;
655 Float_t etaProng=0., ptProng=0.;
657 for (Int_t iProng =0; iProng<fProngs; iProng++){
659 trackCuts[iProng]->GetEtaRange(etaCutMin, etaCutMax);
660 trackCuts[iProng]->GetPtRange(ptCutMin, ptCutMax);
661 etaProng = GetEtaProng(iProng);
662 ptProng = GetPtProng(iProng);
664 Bool_t acceptanceProng = (etaProng>etaCutMin && etaProng<etaCutMax && ptProng>ptCutMin && ptProng<ptCutMax);
665 if (!acceptanceProng) {
666 AliDebug(2,"At least one reconstructed prong isn't in the acceptance\n");
674 //___________________________________________________________
676 Bool_t AliCFVertexingHF::FillUnfoldingMatrix(Double_t fill[4]) const
679 // filling the unfolding matrix
682 if(fmcPartCandidate){
684 fill[0] = GetPtCand();
685 fill[1] = GetYCand();
687 fill[2] = fmcPartCandidate->Pt();
688 fill[3] = fmcPartCandidate->Y();
695 //___________________________________________________________
697 Int_t AliCFVertexingHF::CheckReflexion(Char_t isSign)
700 // check for reflexion (particle/antiparticle)
703 Int_t mcLabel = GetMCLabel();
706 AliDebug(2,"No MC particle found");
710 fmcPartCandidate = (AliAODMCParticle*)fmcArray->At(mcLabel);
711 if (!fmcPartCandidate){
712 AliWarning("Could not find associated MC in AOD MC tree");
717 if(fmcPartCandidate->GetPdgCode()>0) {
718 if (isSign == 1){ // I ask for antiparticle only
719 AliDebug(2,"candidate is particle, I ask for antiparticle only");
722 return 1; // particle
724 else if(fmcPartCandidate->GetPdgCode()<0) {
725 if (isSign == 0){ // I ask for particle only
726 AliDebug(2,"candidate is antiparticle, I ask for particle only");
729 return 2; // antiparticle
731 else return 0; // ....shouldn't be...
734 //___________________________________________________________
736 Bool_t AliCFVertexingHF::SetLabelArray()
739 // setting the label arrays
742 Bool_t bLabelArray = kFALSE;
744 fLabelArray = new Int_t[fProngs];
746 AliAODMCParticle *mcPartDaughter;
747 Int_t label0 = fmcPartCandidate->GetDaughter(0);
748 Int_t label1 = fmcPartCandidate->GetDaughter(1);
749 AliDebug(2,Form("label0 = %d, label1 = %d",label0,label1));
750 if (label1==0 || label0 == 0){
751 AliDebug(2, Form("The MC particle doesn't have correct daughters, skipping!!"));
752 delete [] fLabelArray;
757 if (label1 - label0 == fProngs-1){
758 for (Int_t iProng = 0; iProng<fProngs; iProng++){
759 mcPartDaughter = dynamic_cast<AliAODMCParticle*>(fmcArray->At(label0+iProng));
761 fLabelArray[iProng] = mcPartDaughter->GetLabel();
764 AliError("Failed casting the daughter particle, returning a NULL label array");
765 delete [] fLabelArray;
772 // resonant decay channel
773 else if (label1 - label0 == fProngs-2 && fProngs > 2){
774 Int_t labelFirstDau = fmcPartCandidate->GetDaughter(0);
775 Int_t foundDaughters = 0;
776 for(Int_t iDau=0; iDau<fProngs-1; iDau++){
777 Int_t iLabelDau = labelFirstDau+iDau;
778 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iLabelDau));
780 AliError("Wrong particle type in fmcArray");
781 delete [] fLabelArray;
785 Int_t pdgCode=TMath::Abs(part->GetPdgCode());
786 if(pdgCode==211 || pdgCode==321 || pdgCode==2212){
788 fLabelArray[foundDaughters] = part->GetLabel();
792 AliError("Error while casting particle! returning a NULL array");
793 delete [] fLabelArray;
799 Int_t nDauRes=part->GetNDaughters();
801 delete [] fLabelArray;
805 Int_t labelFirstDauRes = part->GetDaughter(0);
806 for(Int_t iDauRes=0; iDauRes<nDauRes; iDauRes++){
807 Int_t iLabelDauRes = labelFirstDauRes+iDauRes;
808 AliAODMCParticle* dauRes = dynamic_cast<AliAODMCParticle*>(fmcArray->At(iLabelDauRes));
810 fLabelArray[foundDaughters] = dauRes->GetLabel();
814 AliError("Error while casting resonant daughter! returning a NULL array");
815 delete [] fLabelArray;
822 if (foundDaughters != fProngs){
823 delete [] fLabelArray;
828 // wrong correspondance label <--> prongs
830 delete [] fLabelArray;
834 SetAccCut(); // setting the pt and eta acceptance cuts
839 //___________________________________________________________
841 void AliCFVertexingHF::SetPtAccCut(Float_t* ptAccCut)
844 // setting the pt cut to be used in the Acceptance steps (MC+Reco)
848 for (Int_t iP=0; iP<fProngs; iP++){
849 fPtAccCut[iP]=ptAccCut[iP];
857 //___________________________________________________________
859 void AliCFVertexingHF::SetEtaAccCut(Float_t* etaAccCut)
862 // setting the eta cut to be used in the Acceptance steps (MC+Reco)
866 for (Int_t iP=0; iP<fProngs; iP++){
867 fEtaAccCut[iP]=etaAccCut[iP];
872 //___________________________________________________________
874 void AliCFVertexingHF::SetAccCut(Float_t* ptAccCut, Float_t* etaAccCut)
877 // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
881 for (Int_t iP=0; iP<fProngs; iP++){
882 fPtAccCut[iP]=ptAccCut[iP];
883 fEtaAccCut[iP]=etaAccCut[iP];
889 //___________________________________________________________
891 void AliCFVertexingHF::SetAccCut()
894 // setting the pt and eta cut to be used in the Acceptance steps (MC+Reco)
898 for (Int_t iP=0; iP<fProngs; iP++){