[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliVertexingHFUtils.cxx

/**************************************************************************
 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

#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$ */

///////////////////////////////////////////////////////////////////
//                                                               //
// Class with functions useful for different D2H analyses        //
// - event plane resolution                                      //
// - <pt> calculation with side band subtraction                 //
// - tracklet multiplicity calculation                            //
// Origin: F.Prino, Torino, prino@to.infn.it                     //
//                                                               //
///////////////////////////////////////////////////////////////////

ClassImp(AliVertexingHFUtils)

//______________________________________________________________________
AliVertexingHFUtils::AliVertexingHFUtils():TObject(),
  fK(1),
  fSubRes(1.),
  fMinEtaForTracklets(-1.),
  fMaxEtaForTracklets(1.)
{
  // Default contructor
}


//______________________________________________________________________
AliVertexingHFUtils::AliVertexingHFUtils(Int_t k):
  TObject(),
  fK(k),
  fSubRes(1.),
  fMinEtaForTracklets(-1.),
  fMaxEtaForTracklets(1.)
{
  // Standard constructor
}


//______________________________________________________________________
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
  Double_t c[5];
  if(k==1){ 
    c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
  }
  else if(k==2){
    c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
  }
  else return -1;
  return c[0]*x+c[1]*x*x+c[2]*x*x*x+c[3]*x*x*x*x+c[4]*x*x*x*x*x;
}

//______________________________________________________________________
Double_t AliVertexingHFUtils:: ResolK1(Double_t x){
  return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4));
}


//______________________________________________________________________
Double_t AliVertexingHFUtils::FindChi(Double_t res,  Int_t k){
  // compute chi variable (=v2*sqrt(N)) from external values

  Double_t x1=0;
  Double_t x2=15;
  Double_t y1,y2;
  if(k==1){
    y1=ResolK1(x1)-res;
    y2=ResolK1(x2)-res;
  }
  else if(k==2){
    y1=Pol(x1,2)-res;
    y2=Pol(x2,2)-res;
  }
  else return -1;

  if(y1*y2>0) return -1;
  if(y1==0) return y1;
  if(y2==0) return y2;
  Double_t xmed,ymed;
  Int_t jiter=0;
  while((x2-x1)>0.0001){
    xmed=0.5*(x1+x2);
    if(k==1){
      y1=ResolK1(x1)-res;
      y2=ResolK1(x2)-res;
      ymed=ResolK1(xmed)-res;
    }
    else if(k==2){
      y1=Pol(x1,2)-res;
      y2=Pol(x2,2)-res;
      ymed=Pol(xmed,2)-res;
    }
    else return -1;
    if((y1*ymed)<0) x2=xmed;
    if((y2*ymed)<0)x1=xmed;
    if(ymed==0) return xmed;
    jiter++;
  }
  return 0.5*(x1+x2);
}

//______________________________________________________________________
Double_t AliVertexingHFUtils::GetFullEvResol(Double_t resSub, Int_t k){
  // computes event plane resolution starting from sub event resolution
  Double_t chisub=FindChi(resSub,k);
  Double_t chifull=chisub*TMath::Sqrt(2);
  if(k==1) return ResolK1(chifull);
  else if(k==2) return Pol(chifull,2);
  else{
    printf("k should be <=2\n");
    return 1.;
  }
}

//______________________________________________________________________
Double_t AliVertexingHFUtils::GetFullEvResol(const TH1F* hSubEvCorr, Int_t k){
  // computes event plane resolution starting from sub event correlation histogram
  if(!hSubEvCorr) return 1.;
  Double_t resSub=GetSubEvResol(hSubEvCorr);
  return GetFullEvResol(resSub,k);
}
//______________________________________________________________________
Double_t AliVertexingHFUtils::GetFullEvResolLowLim(const TH1F* hSubEvCorr, Int_t k){
  // computes low limit event plane resolution starting from sub event correlation histogram
  if(!hSubEvCorr) return 1.;
  Double_t resSub=GetSubEvResolLowLim(hSubEvCorr);
  printf("%f\n",resSub);
  return GetFullEvResol(resSub,k);  
}
//______________________________________________________________________
Double_t AliVertexingHFUtils::GetFullEvResolHighLim(const TH1F* hSubEvCorr, Int_t k){
  // computes high limit event plane resolution starting from sub event correlation histogram
  if(!hSubEvCorr) return 1.;
  Double_t resSub=GetSubEvResolHighLim(hSubEvCorr);
  printf("%f\n",resSub);
  return GetFullEvResol(resSub,k);  
}
//______________________________________________________________________
Int_t AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(AliAODEvent* ev, Double_t mineta, Double_t maxeta){
  // counts tracklets in given eta range
  AliAODTracklets* tracklets=ev->GetTracklets();
  Int_t nTr=tracklets->GetNumberOfTracklets();
  Int_t count=0;
  for(Int_t iTr=0; iTr<nTr; iTr++){
    Double_t theta=tracklets->GetTheta(iTr);
    Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
    if(eta>mineta && eta<maxeta) count++;
  }
  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

  //Make 2D histos in the desired pt range
  Int_t start=hMassD->FindBin(ptmin);
  Int_t end=hMassD->FindBin(ptmax)-1;
  const Int_t nx=end-start;
  TH2F *hMassDpt=new TH2F("hptmass","hptmass",nx,ptmin,ptmax,hMassD->GetNbinsY(),hMassD->GetYaxis()->GetBinLowEdge(1),hMassD->GetYaxis()->GetBinLowEdge(hMassD->GetNbinsY())+hMassD->GetYaxis()->GetBinWidth(hMassD->GetNbinsY()));
  for(Int_t ix=start;ix<end;ix++){
    for(Int_t iy=1;iy<=hMassD->GetNbinsY();iy++){
      hMassDpt->SetBinContent(ix-start+1,iy,hMassD->GetBinContent(ix,iy));
      hMassDpt->SetBinError(ix-start+1,iy,hMassD->GetBinError(ix,iy));
    }
  }

  Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit;
  Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit;
  Int_t minBinSig=hMassD->GetYaxis()->FindBin(minMassSig);
  Int_t maxBinSig=hMassD->GetYaxis()->FindBin(maxMassSig);
  Double_t minMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(minBinSig);
  Double_t maxMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinSig)+hMassD->GetYaxis()->GetBinWidth(maxBinSig);
  //  printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);

  Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
  Int_t minBinBkgLow=2;
  Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow);
  Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow);
  Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow);
  Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
  Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi);
  Int_t maxBinBkgHi=hMassD->GetNbinsY()-1;
  Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi);
  Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi);
  //  printf("BKG Fit Limits = %f %f  && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);

  Double_t bkgSig=funcB2->Integral(minMassSigBin,maxMassSigBin);
  Double_t bkgLow=funcB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
  Double_t bkgHi=funcB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
  //  printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);

  TH1F* hMptBkgLo=(TH1F*)hMassDpt->ProjectionX("hPtBkgLoBin",minBinBkgLow,maxBinBkgLow);
  TH1F* hMptBkgHi=(TH1F*)hMassDpt->ProjectionX("hPtBkgHiBin",minBinBkgHi,maxBinBkgHi);
  TH1F* hMptSigReg=(TH1F*)hMassDpt->ProjectionX("hCPtBkgSigBin",minBinSig,maxBinSig);

  hMptBkgLo->Rebin(rebin);
  hMptBkgHi->Rebin(rebin);
  hMptSigReg->Rebin(rebin);

  hMptBkgLo->Sumw2();
  hMptBkgHi->Sumw2();
  TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin");
  hMptBkgLoScal->Scale(bkgSig/bkgLow);
  TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin");
  hMptBkgHiScal->Scale(bkgSig/bkgHi);

  TH1F* hMptBkgAver=0x0;
  hMptBkgAver=(TH1F*)hMptBkgLoScal->Clone("hPtBkgAverBin");
  hMptBkgAver->Add(hMptBkgHiScal);
  hMptBkgAver->Scale(0.5);
  TH1F* hMptSig=(TH1F*)hMptSigReg->Clone("hCPtSigBin");
  hMptSig->Add(hMptBkgAver,-1.);   
 
  averagePt = hMptSig->GetMean();
  errorPt = hMptSig->GetMeanError();

  delete hMptBkgLo;
  delete hMptBkgHi;
  delete hMptSigReg;
  delete hMptBkgLoScal;
  delete hMptBkgHiScal;
  delete hMptBkgAver;
  delete hMassDpt;
  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;
 }
Commit	Line	Data
a6c5a2e9	1	/**************************************************************************
	2	* Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	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	**************************************************************************/
	15
	16	#include <TMath.h>
686e4710	17	#include <TRandom.h>
b176db3b	18	#include <TProfile.h>
	19	#include <TClonesArray.h>
	20	#include <TH1F.h>
	21	#include <TH2F.h>
	22	#include <TF1.h>
a6c5a2e9	23	#include "AliAODEvent.h"
b176db3b	24	#include "AliAODMCHeader.h"
7f51a627	25	#include "AliGenEventHeader.h"
25504150	26	#include "AliAODMCParticle.h"
b176db3b	27	#include "AliAODRecoDecayHF.h"
a6c5a2e9	28	#include "AliVertexingHFUtils.h"
a6c5a2e9	29
d947ea9c	30	/* $Id$ */
a6c5a2e9	31
	32	///////////////////////////////////////////////////////////////////
	33	// //
	34	// Class with functions useful for different D2H analyses //
	35	// - event plane resolution //
	36	// - <pt> calculation with side band subtraction //
25504150	37	// - tracklet multiplicity calculation //
a6c5a2e9	38	// Origin: F.Prino, Torino, prino@to.infn.it //
	39	// //
	40	///////////////////////////////////////////////////////////////////
	41
	42	ClassImp(AliVertexingHFUtils)
	43
	44	//______________________________________________________________________
	45	AliVertexingHFUtils::AliVertexingHFUtils():TObject(),
	46	fK(1),
	47	fSubRes(1.),
	48	fMinEtaForTracklets(-1.),
	49	fMaxEtaForTracklets(1.)
	50	{
	51	// Default contructor
	52	}
	53
	54
	55	//______________________________________________________________________
	56	AliVertexingHFUtils::AliVertexingHFUtils(Int_t k):
	57	TObject(),
	58	fK(k),
	59	fSubRes(1.),
	60	fMinEtaForTracklets(-1.),
	61	fMaxEtaForTracklets(1.)
	62	{
	63	// Standard constructor
	64	}
	65
	66
5d1396f9	67	//______________________________________________________________________
	68	void AliVertexingHFUtils::ComputeSignificance(Double_t signal, Double_t errsignal, Double_t background, Double_t errbackground, Double_t &significance,Double_t &errsignificance){
	69	// calculate significance from S, B and errors
	70
	71
	72	Double_t errSigSq=errsignal*errsignal;
	73	Double_t errBkgSq=errbackground*errbackground;
	74	Double_t sigPlusBkg=signal+background;
	75	if(sigPlusBkg>0. && signal>0.){
	76	significance = signal/TMath::Sqrt(signal+background);
	77	errsignificance = significanceTMath::Sqrt((errSigSq+errBkgSq)/(4.sigPlusBkgsigPlusBkg)+(background/sigPlusBkg)errSigSq/signal/signal);
	78	}else{
	79	significance=0.;
	80	errsignificance=0.;
	81	}
	82	return;
	83
	84	}
a6c5a2e9	85	//______________________________________________________________________
	86	Double_t AliVertexingHFUtils::Pol(Double_t x, Int_t k){
	87	// compute chi from polynomial approximation
	88	Double_t c[5];
	89	if(k==1){
	90	c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
	91	}
	92	else if(k==2){
	93	c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
	94	}
	95	else return -1;
	96	return c[0]x+c[1]xx+c[2]xxx+c[3]xxxx+c[4]xxxx*x;
	97	}
	98
	99	//______________________________________________________________________
	100	Double_t AliVertexingHFUtils:: ResolK1(Double_t x){
	101	return TMath::Sqrt(TMath::Pi()/8)xTMath::Exp(-xx/4)(TMath::BesselI0(xx/4)+TMath::BesselI1(xx/4));
	102	}
	103
	104
	105	//______________________________________________________________________
	106	Double_t AliVertexingHFUtils::FindChi(Double_t res, Int_t k){
	107	// compute chi variable (=v2*sqrt(N)) from external values
	108
	109	Double_t x1=0;
	110	Double_t x2=15;
	111	Double_t y1,y2;
	112	if(k==1){
	113	y1=ResolK1(x1)-res;
	114	y2=ResolK1(x2)-res;
	115	}
	116	else if(k==2){
	117	y1=Pol(x1,2)-res;
	118	y2=Pol(x2,2)-res;
	119	}
	120	else return -1;
	121
	122	if(y1*y2>0) return -1;
	123	if(y1==0) return y1;
	124	if(y2==0) return y2;
	125	Double_t xmed,ymed;
	126	Int_t jiter=0;
	127	while((x2-x1)>0.0001){
	128	xmed=0.5*(x1+x2);
	129	if(k==1){
	130	y1=ResolK1(x1)-res;
	131	y2=ResolK1(x2)-res;
	132	ymed=ResolK1(xmed)-res;
	133	}
	134	else if(k==2){
	135	y1=Pol(x1,2)-res;
	136	y2=Pol(x2,2)-res;
	137	ymed=Pol(xmed,2)-res;
	138	}
	139	else return -1;
	140	if((y1*ymed)<0) x2=xmed;
	141	if((y2*ymed)<0)x1=xmed;
	142	if(ymed==0) return xmed;
	143	jiter++;
	144	}
	145	return 0.5*(x1+x2);
	146	}
	147
	148	//______________________________________________________________________
149	Double_t AliVertexingHFUtils::GetFullEvResol(Double_t resSub, Int_t k){
150	// computes event plane resolution starting from sub event resolution
151	Double_t chisub=FindChi(resSub,k);
152	Double_t chifull=chisub*TMath::Sqrt(2);
153	if(k==1) return ResolK1(chifull);
154	else if(k==2) return Pol(chifull,2);
155	else{
156	printf("k should be <=2\n");
157	return 1.;
158	}
159	}
160
161	//______________________________________________________________________
162	Double_t AliVertexingHFUtils::GetFullEvResol(const TH1F* hSubEvCorr, Int_t k){
163	// computes event plane resolution starting from sub event correlation histogram
164	if(!hSubEvCorr) return 1.;
165	Double_t resSub=GetSubEvResol(hSubEvCorr);
166	return GetFullEvResol(resSub,k);
167	}
168	//______________________________________________________________________
169	Double_t AliVertexingHFUtils::GetFullEvResolLowLim(const TH1F* hSubEvCorr, Int_t k){
170	// computes low limit event plane resolution starting from sub event correlation histogram
171	if(!hSubEvCorr) return 1.;
172	Double_t resSub=GetSubEvResolLowLim(hSubEvCorr);
173	printf("%f\n",resSub);
174	return GetFullEvResol(resSub,k);
175	}
176	//______________________________________________________________________
177	Double_t AliVertexingHFUtils::GetFullEvResolHighLim(const TH1F* hSubEvCorr, Int_t k){
178	// computes high limit event plane resolution starting from sub event correlation histogram
179	if(!hSubEvCorr) return 1.;
180	Double_t resSub=GetSubEvResolHighLim(hSubEvCorr);
181	printf("%f\n",resSub);
182	return GetFullEvResol(resSub,k);
183	}
184	//______________________________________________________________________
185	Int_t AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(AliAODEvent* ev, Double_t mineta, Double_t maxeta){
186	// counts tracklets in given eta range
187	AliAODTracklets* tracklets=ev->GetTracklets();
188	Int_t nTr=tracklets->GetNumberOfTracklets();
189	Int_t count=0;
190	for(Int_t iTr=0; iTr<nTr; iTr++){
191	Double_t theta=tracklets->GetTheta(iTr);
192	Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
193	if(eta>mineta && eta<maxeta) count++;
194	}
195	return count;
196	}
86d84fa4	197	//______________________________________________________________________
25504150	198	Int_t AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
	199	// counts generated particles in fgiven eta range
	200
	201	Int_t nChargedMC=0;
	202	for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
	203	AliAODMCParticle part=(AliAODMCParticle)arrayMC->UncheckedAt(i);
	204	Int_t charge = part->Charge();
	205	Double_t eta = part->Eta();
	206	if(charge!=0 && eta>mineta && eta<maxeta) nChargedMC++;
	207	}
	208	return nChargedMC;
	209	}
	210	//______________________________________________________________________
	211	Int_t AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
	212	// counts generated primary particles in given eta range
	213
	214	Int_t nChargedMC=0;
	215	for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
	216	AliAODMCParticle part=(AliAODMCParticle)arrayMC->UncheckedAt(i);
	217	Int_t charge = part->Charge();
	218	Double_t eta = part->Eta();
	219	if(charge!=0 && eta>mineta && eta<maxeta){
	220	if(part->IsPrimary())nChargedMC++;
	221	}
	222	}
	223	return nChargedMC;
	224	}
	225	//______________________________________________________________________
	226	Int_t AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
	227	// counts generated primary particles in given eta range
	228
	229	Int_t nChargedMC=0;
	230	for(Int_t i=0;i<arrayMC->GetEntriesFast();i++){
	231	AliAODMCParticle part=(AliAODMCParticle)arrayMC->UncheckedAt(i);
	232	Int_t charge = part->Charge();
	233	Double_t eta = part->Eta();
	234	if(charge!=0 && eta>mineta && eta<maxeta){
	235	if(part->IsPhysicalPrimary())nChargedMC++;
	236	}
	237	}
	238	return nChargedMC;
	239	}
	240	//______________________________________________________________________
86d84fa4	241	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){
	242
	243	// Compute <pt> from 2D histogram M vs pt
	244
	245	//Make 2D histos in the desired pt range
	246	Int_t start=hMassD->FindBin(ptmin);
	247	Int_t end=hMassD->FindBin(ptmax)-1;
	248	const Int_t nx=end-start;
	249	TH2F *hMassDpt=new TH2F("hptmass","hptmass",nx,ptmin,ptmax,hMassD->GetNbinsY(),hMassD->GetYaxis()->GetBinLowEdge(1),hMassD->GetYaxis()->GetBinLowEdge(hMassD->GetNbinsY())+hMassD->GetYaxis()->GetBinWidth(hMassD->GetNbinsY()));
	250	for(Int_t ix=start;ix<end;ix++){
	251	for(Int_t iy=1;iy<=hMassD->GetNbinsY();iy++){
	252	hMassDpt->SetBinContent(ix-start+1,iy,hMassD->GetBinContent(ix,iy));
	253	hMassDpt->SetBinError(ix-start+1,iy,hMassD->GetBinError(ix,iy));
	254	}
	255	}
	256
	257	Double_t minMassSig=massFromFit-sigmaRangeForSig*sigmaFromFit;
	258	Double_t maxMassSig=massFromFit+sigmaRangeForSig*sigmaFromFit;
	259	Int_t minBinSig=hMassD->GetYaxis()->FindBin(minMassSig);
	260	Int_t maxBinSig=hMassD->GetYaxis()->FindBin(maxMassSig);
	261	Double_t minMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(minBinSig);
	262	Double_t maxMassSigBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinSig)+hMassD->GetYaxis()->GetBinWidth(maxBinSig);
	263	// printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
	264
	265	Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
	266	Int_t minBinBkgLow=2;
	267	Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow);
	268	Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow);
	269	Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow);
	270	Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
	271	Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi);
	272	Int_t maxBinBkgHi=hMassD->GetNbinsY()-1;
	273	Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi);
	274	Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi);
	275	// printf("BKG Fit Limits = %f %f && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
	276
	277	Double_t bkgSig=funcB2->Integral(minMassSigBin,maxMassSigBin);
	278	Double_t bkgLow=funcB2->Integral(minMassBkgLowBin,maxMassBkgLowBin);
	279	Double_t bkgHi=funcB2->Integral(minMassBkgHiBin,maxMassBkgHiBin);
	280	// printf("Background integrals = %f %f %f\n",bkgLow,bkgSig,bkgHi);
	281
	282	TH1F* hMptBkgLo=(TH1F*)hMassDpt->ProjectionX("hPtBkgLoBin",minBinBkgLow,maxBinBkgLow);
	283	TH1F* hMptBkgHi=(TH1F*)hMassDpt->ProjectionX("hPtBkgHiBin",minBinBkgHi,maxBinBkgHi);
	284	TH1F* hMptSigReg=(TH1F*)hMassDpt->ProjectionX("hCPtBkgSigBin",minBinSig,maxBinSig);
	285
	286	hMptBkgLo->Rebin(rebin);
	287	hMptBkgHi->Rebin(rebin);
	288	hMptSigReg->Rebin(rebin);
	289
	290	hMptBkgLo->Sumw2();
	291	hMptBkgHi->Sumw2();
	292	TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin");
	293	hMptBkgLoScal->Scale(bkgSig/bkgLow);
	294	TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin");
	295	hMptBkgHiScal->Scale(bkgSig/bkgHi);
	296
	297	TH1F* hMptBkgAver=0x0;
	298	hMptBkgAver=(TH1F*)hMptBkgLoScal->Clone("hPtBkgAverBin");
	299	hMptBkgAver->Add(hMptBkgHiScal);
	300	hMptBkgAver->Scale(0.5);
	301	TH1F* hMptSig=(TH1F*)hMptSigReg->Clone("hCPtSigBin");
	302	hMptSig->Add(hMptBkgAver,-1.);
	303
	304	averagePt = hMptSig->GetMean();
305	errorPt = hMptSig->GetMeanError();
306
307	delete hMptBkgLo;
308	delete hMptBkgHi;
309	delete hMptSigReg;
310	delete hMptBkgLoScal;
311	delete hMptBkgHiScal;
312	delete hMptBkgAver;
313	delete hMassDpt;
314	delete hMptSig;
315
15917f37	316	}
	317	//____________________________________________________________________________
	318	Double_t AliVertexingHFUtils::GetTrueImpactParameterDzero(AliAODMCHeader mcHeader, TClonesArray arrayMC, AliAODMCParticle *partD) {
	319	// true impact parameter calculation for Dzero
	320
	321	if(!partD \|\| !arrayMC \|\| !mcHeader) return 99999.;
	322	Int_t code=partD->GetPdgCode();
	323	if(TMath::Abs(code)!=421) return 99999.;
	324
	325	Double_t vtxTrue[3];
	326	mcHeader->GetVertex(vtxTrue);
	327	Double_t origD[3];
	328	partD->XvYvZv(origD);
	329	Short_t charge=partD->Charge();
	330	Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
	331	for(Int_t iDau=0; iDau<2; iDau++){
	332	pXdauTrue[iDau]=0.;
	333	pYdauTrue[iDau]=0.;
	334	pZdauTrue[iDau]=0.;
	335	}
	336
	337	Int_t nDau=partD->GetNDaughters();
	338	Int_t labelFirstDau = partD->GetDaughter(0);
	339	if(nDau==2){
	340	for(Int_t iDau=0; iDau<2; iDau++){
	341	Int_t ind = labelFirstDau+iDau;
	342	AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
	343	if(!part){
	344	printf("Daughter particle not found in MC array");
	345	return 99999.;
	346	}
	347	pXdauTrue[iDau]=part->Px();
	348	pYdauTrue[iDau]=part->Py();
	349	pZdauTrue[iDau]=part->Pz();
	350	}
	351	}else{
	352	return 99999.;
	353	}
	354
48efc6ad	355	Double_t d0dummy[2]={0.,0.};
15917f37	356	AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
	357	return aodDvsMC.ImpParXY();
	358
86d84fa4	359	}
b176db3b	360	//____________________________________________________________________________
	361	Double_t AliVertexingHFUtils::GetTrueImpactParameterDplus(AliAODMCHeader mcHeader, TClonesArray arrayMC, AliAODMCParticle *partD) {
	362	// true impact parameter calculation for Dplus
	363
	364	if(!partD \|\| !arrayMC \|\| !mcHeader) return 99999.;
	365	Int_t code=partD->GetPdgCode();
	366	if(TMath::Abs(code)!=411) return 99999.;
	367
	368	Double_t vtxTrue[3];
	369	mcHeader->GetVertex(vtxTrue);
	370	Double_t origD[3];
	371	partD->XvYvZv(origD);
	372	Short_t charge=partD->Charge();
	373	Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
	374	for(Int_t iDau=0; iDau<3; iDau++){
	375	pXdauTrue[iDau]=0.;
	376	pYdauTrue[iDau]=0.;
	377	pZdauTrue[iDau]=0.;
	378	}
	379
	380	Int_t nDau=partD->GetNDaughters();
	381	Int_t labelFirstDau = partD->GetDaughter(0);
	382	if(nDau==3){
	383	for(Int_t iDau=0; iDau<3; iDau++){
	384	Int_t ind = labelFirstDau+iDau;
	385	AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
	386	if(!part){
	387	printf("Daughter particle not found in MC array");
	388	return 99999.;
	389	}
	390	pXdauTrue[iDau]=part->Px();
	391	pYdauTrue[iDau]=part->Py();
	392	pZdauTrue[iDau]=part->Pz();
	393	}
	394	}else if(nDau==2){
	395	Int_t theDau=0;
	396	for(Int_t iDau=0; iDau<2; iDau++){
	397	Int_t ind = labelFirstDau+iDau;
	398	AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
	399	if(!part){
	400	printf("Daughter particle not found in MC array");
	401	return 99999.;
	402	}
	403	Int_t pdgCode=TMath::Abs(part->GetPdgCode());
	404	if(pdgCode==211 \|\| pdgCode==321){
	405	pXdauTrue[theDau]=part->Px();
	406	pYdauTrue[theDau]=part->Py();
	407	pZdauTrue[theDau]=part->Pz();
	408	++theDau;
	409	}else{
	410	Int_t nDauRes=part->GetNDaughters();
	411	if(nDauRes==2){
	412	Int_t labelFirstDauRes = part->GetDaughter(0);
	413	for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
	414	Int_t indDR = labelFirstDauRes+iDauRes;
	415	AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
	416	if(!partDR){
	417	printf("Daughter particle not found in MC array");
	418	return 99999.;
	419	}
	420
	421	Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
	422	if(pdgCodeDR==211 \|\| pdgCodeDR==321){
	423	pXdauTrue[theDau]=partDR->Px();
424	pYdauTrue[theDau]=partDR->Py();
425	pZdauTrue[theDau]=partDR->Pz();
426	++theDau;
427	}
428	}
429	}
430	}
431	}
432	if(theDau!=3){
433	printf("Wrong number of decay prongs");
434	return 99999.;
435	}
436	}
437
438	Double_t d0dummy[3]={0.,0.,0.};
439	AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
440	return aodDvsMC.ImpParXY();
441
442	}
443
444
445
686e4710	446	//____________________________________________________________________________
	447	Double_t AliVertexingHFUtils::GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult) {
	448	//
	449	// Correct the number of accepted tracklets based on a TProfile Hist
	450	//
	451	//
	452
	453	if(TMath::Abs(vtxZ)>10.0){
	454	printf("ERROR: Z vertex out of range for correction of multiplicity\n");
	455	return uncorrectedNacc;
	456	}
	457
	458	if(!estimatorAvg){
	459	printf("ERROR: Missing TProfile for correction of multiplicity\n");
	460	return uncorrectedNacc;
	461	}
	462
	463	Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));
	464
	465	Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);
	466
	467	Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM));
	468
	469	return correctedNacc;
	470	}
7f51a627	471	//______________________________________________________________________
	472	TString AliVertexingHFUtils::GetGenerator(Int_t label, AliAODMCHeader* header){
	473	Int_t nsumpart=0;
	474	TList *lh=header->GetCocktailHeaders();
	475	Int_t nh=lh->GetEntries();
	476	for(Int_t i=0;i<nh;i++){
	477	AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
	478	TString genname=gh->GetName();
	479	Int_t npart=gh->NProduced();
	480	if(label>=nsumpart && label<(nsumpart+npart)) return genname;
	481	nsumpart+=npart;
	482	}
	483	TString empty="";
	484	return empty;
	485	}