**************************************************************************/
#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 "AliGenEventHeader.h"
+#include "AliAODMCParticle.h"
+#include "AliAODRecoDecayHF.h"
#include "AliVertexingHFUtils.h"
-/* $Id: $ */
+/* $Id$ */
///////////////////////////////////////////////////////////////////
// //
// Class with functions useful for different D2H analyses //
// - event plane resolution //
// - <pt> calculation with side band subtraction //
-// - tracklet multiplcity calculation //
+// - tracklet multiplicity calculation //
// Origin: F.Prino, Torino, prino@to.infn.it //
// //
///////////////////////////////////////////////////////////////////
}
+//______________________________________________________________________
+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
return count;
}
//______________________________________________________________________
+Int_t AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
+ // counts generated particles in fgiven eta range
+
+ Int_t nChargedMC=0;
+ for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
+ AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
+ Int_t charge = part->Charge();
+ Double_t eta = part->Eta();
+ if(charge!=0 && eta>mineta && eta<maxeta) nChargedMC++;
+ }
+ return nChargedMC;
+}
+//______________________________________________________________________
+Int_t AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
+ // counts generated primary particles in given eta range
+
+ Int_t nChargedMC=0;
+ for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
+ AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
+ Int_t charge = part->Charge();
+ Double_t eta = part->Eta();
+ if(charge!=0 && eta>mineta && eta<maxeta){
+ if(part->IsPrimary())nChargedMC++;
+ }
+ }
+ return nChargedMC;
+}
+//______________________________________________________________________
+Int_t AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
+ // counts generated primary particles in given eta range
+
+ Int_t nChargedMC=0;
+ for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
+ AliAODMCParticle *part=(AliAODMCParticle*)arrayMC->UncheckedAt(i);
+ Int_t charge = part->Charge();
+ Double_t eta = part->Eta();
+ if(charge!=0 && eta>mineta && eta<maxeta){
+ if(part->IsPhysicalPrimary())nChargedMC++;
+ }
+ }
+ return nChargedMC;
+}
+//______________________________________________________________________
void AliVertexingHFUtils::AveragePt(Float_t& averagePt, Float_t& errorPt,Float_t ptmin,Float_t ptmax, TH2F* hMassD, Float_t massFromFit, Float_t sigmaFromFit, TF1* funcB2, Float_t sigmaRangeForSig,Float_t sigmaRangeForBkg, Int_t rebin){
// Compute <pt> from 2D histogram M vs pt
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;
+}
+//______________________________________________________________________
+TString AliVertexingHFUtils::GetGenerator(Int_t label, AliAODMCHeader* header){
+ Int_t nsumpart=0;
+ TList *lh=header->GetCocktailHeaders();
+ Int_t nh=lh->GetEntries();
+ for(Int_t i=0;i<nh;i++){
+ AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
+ TString genname=gh->GetName();
+ Int_t npart=gh->NProduced();
+ if(label>=nsumpart && label<(nsumpart+npart)) return genname;
+ nsumpart+=npart;
+ }
+ TString empty="";
+ return empty;
+ }
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff