**************************************************************************/
#include <TMath.h>
+#include <TRandom.h>
+#include <TProfile.h>
+#include <TClonesArray.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TF1.h>
#include "AliAODEvent.h"
+#include "AliAODMCHeader.h"
#include "AliAODMCParticle.h"
+#include "AliAODRecoDecayHF.h"
#include "AliVertexingHFUtils.h"
-/* $Id: $ */
+/* $Id$ */
///////////////////////////////////////////////////////////////////
// //
}
+//______________________________________________________________________
+void AliVertexingHFUtils::ComputeSignificance(Double_t signal, Double_t errsignal, Double_t background, Double_t errbackground, Double_t &significance,Double_t &errsignificance){
+ // calculate significance from S, B and errors
+
+
+ Double_t errSigSq=errsignal*errsignal;
+ Double_t errBkgSq=errbackground*errbackground;
+ Double_t sigPlusBkg=signal+background;
+ if(sigPlusBkg>0. && signal>0.){
+ significance = signal/TMath::Sqrt(signal+background);
+ errsignificance = significance*TMath::Sqrt((errSigSq+errBkgSq)/(4.*sigPlusBkg*sigPlusBkg)+(background/sigPlusBkg)*errSigSq/signal/signal);
+ }else{
+ significance=0.;
+ errsignificance=0.;
+ }
+ return;
+
+}
//______________________________________________________________________
Double_t AliVertexingHFUtils::Pol(Double_t x, Int_t k){
// compute chi from polynomial approximation
delete hMptSig;
}
+//____________________________________________________________________________
+Double_t AliVertexingHFUtils::GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
+ // true impact parameter calculation for Dzero
+
+ if(!partD || !arrayMC || !mcHeader) return 99999.;
+ Int_t code=partD->GetPdgCode();
+ if(TMath::Abs(code)!=421) return 99999.;
+
+ Double_t vtxTrue[3];
+ mcHeader->GetVertex(vtxTrue);
+ Double_t origD[3];
+ partD->XvYvZv(origD);
+ Short_t charge=partD->Charge();
+ Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
+ for(Int_t iDau=0; iDau<2; iDau++){
+ pXdauTrue[iDau]=0.;
+ pYdauTrue[iDau]=0.;
+ pZdauTrue[iDau]=0.;
+ }
+
+ Int_t nDau=partD->GetNDaughters();
+ Int_t labelFirstDau = partD->GetDaughter(0);
+ if(nDau==2){
+ for(Int_t iDau=0; iDau<2; iDau++){
+ Int_t ind = labelFirstDau+iDau;
+ AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
+ if(!part){
+ printf("Daughter particle not found in MC array");
+ return 99999.;
+ }
+ pXdauTrue[iDau]=part->Px();
+ pYdauTrue[iDau]=part->Py();
+ pZdauTrue[iDau]=part->Pz();
+ }
+ }else{
+ return 99999.;
+ }
+
+ Double_t d0dummy[2]={0.,0.};
+ AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
+ return aodDvsMC.ImpParXY();
+
+}
+//____________________________________________________________________________
+Double_t AliVertexingHFUtils::GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
+ // true impact parameter calculation for Dplus
+
+ if(!partD || !arrayMC || !mcHeader) return 99999.;
+ Int_t code=partD->GetPdgCode();
+ if(TMath::Abs(code)!=411) return 99999.;
+
+ Double_t vtxTrue[3];
+ mcHeader->GetVertex(vtxTrue);
+ Double_t origD[3];
+ partD->XvYvZv(origD);
+ Short_t charge=partD->Charge();
+ Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
+ for(Int_t iDau=0; iDau<3; iDau++){
+ pXdauTrue[iDau]=0.;
+ pYdauTrue[iDau]=0.;
+ pZdauTrue[iDau]=0.;
+ }
+
+ Int_t nDau=partD->GetNDaughters();
+ Int_t labelFirstDau = partD->GetDaughter(0);
+ if(nDau==3){
+ for(Int_t iDau=0; iDau<3; iDau++){
+ Int_t ind = labelFirstDau+iDau;
+ AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
+ if(!part){
+ printf("Daughter particle not found in MC array");
+ return 99999.;
+ }
+ pXdauTrue[iDau]=part->Px();
+ pYdauTrue[iDau]=part->Py();
+ pZdauTrue[iDau]=part->Pz();
+ }
+ }else if(nDau==2){
+ Int_t theDau=0;
+ for(Int_t iDau=0; iDau<2; iDau++){
+ Int_t ind = labelFirstDau+iDau;
+ AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
+ if(!part){
+ printf("Daughter particle not found in MC array");
+ return 99999.;
+ }
+ Int_t pdgCode=TMath::Abs(part->GetPdgCode());
+ if(pdgCode==211 || pdgCode==321){
+ pXdauTrue[theDau]=part->Px();
+ pYdauTrue[theDau]=part->Py();
+ pZdauTrue[theDau]=part->Pz();
+ ++theDau;
+ }else{
+ Int_t nDauRes=part->GetNDaughters();
+ if(nDauRes==2){
+ Int_t labelFirstDauRes = part->GetDaughter(0);
+ for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
+ Int_t indDR = labelFirstDauRes+iDauRes;
+ AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
+ if(!partDR){
+ printf("Daughter particle not found in MC array");
+ return 99999.;
+ }
+
+ Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
+ if(pdgCodeDR==211 || pdgCodeDR==321){
+ pXdauTrue[theDau]=partDR->Px();
+ pYdauTrue[theDau]=partDR->Py();
+ pZdauTrue[theDau]=partDR->Pz();
+ ++theDau;
+ }
+ }
+ }
+ }
+ }
+ if(theDau!=3){
+ printf("Wrong number of decay prongs");
+ return 99999.;
+ }
+ }
+
+ Double_t d0dummy[3]={0.,0.,0.};
+ AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
+ return aodDvsMC.ImpParXY();
+
+}
+
+
+
+//____________________________________________________________________________
+Double_t AliVertexingHFUtils::GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult) {
+ //
+ // Correct the number of accepted tracklets based on a TProfile Hist
+ //
+ //
+
+ if(TMath::Abs(vtxZ)>10.0){
+ printf("ERROR: Z vertex out of range for correction of multiplicity\n");
+ return uncorrectedNacc;
+ }
+
+ if(!estimatorAvg){
+ printf("ERROR: Missing TProfile for correction of multiplicity\n");
+ return uncorrectedNacc;
+ }
+
+ Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));
+
+ Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);
+
+ Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM));
+
+ return correctedNacc;
+}
+