1 /**************************************************************************
2 * Copyright(c) 2007-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 **************************************************************************/
19 #include <TClonesArray.h>
23 #include "AliAODEvent.h"
24 #include "AliAODMCHeader.h"
25 #include "AliGenEventHeader.h"
26 #include "AliAODMCParticle.h"
27 #include "AliAODRecoDecayHF.h"
28 #include "AliVertexingHFUtils.h"
32 ///////////////////////////////////////////////////////////////////
34 // Class with functions useful for different D2H analyses //
35 // - event plane resolution //
36 // - <pt> calculation with side band subtraction //
37 // - tracklet multiplicity calculation //
38 // Origin: F.Prino, Torino, prino@to.infn.it //
40 ///////////////////////////////////////////////////////////////////
42 ClassImp(AliVertexingHFUtils)
44 //______________________________________________________________________
45 AliVertexingHFUtils::AliVertexingHFUtils():TObject(),
48 fMinEtaForTracklets(-1.),
49 fMaxEtaForTracklets(1.)
55 //______________________________________________________________________
56 AliVertexingHFUtils::AliVertexingHFUtils(Int_t k):
60 fMinEtaForTracklets(-1.),
61 fMaxEtaForTracklets(1.)
63 // Standard constructor
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
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 = significance*TMath::Sqrt((errSigSq+errBkgSq)/(4.*sigPlusBkg*sigPlusBkg)+(background/sigPlusBkg)*errSigSq/signal/signal);
85 //______________________________________________________________________
86 Double_t AliVertexingHFUtils::Pol(Double_t x, Int_t k){
87 // compute chi from polynomial approximation
90 c[0]=0.626657; c[1]=0.; c[2]=-0.09694; c[3]=0.02754; c[4]=-0.002283;
93 c[0]=0.; c[1]=0.25; c[2]=-0.011414; c[3]=-0.034726; c[4]=0.006815;
96 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;
99 //______________________________________________________________________
100 Double_t AliVertexingHFUtils:: ResolK1(Double_t x){
101 return TMath::Sqrt(TMath::Pi()/8)*x*TMath::Exp(-x*x/4)*(TMath::BesselI0(x*x/4)+TMath::BesselI1(x*x/4));
105 //______________________________________________________________________
106 Double_t AliVertexingHFUtils::FindChi(Double_t res, Int_t k){
107 // compute chi variable (=v2*sqrt(N)) from external values
122 if(y1*y2>0) return -1;
127 while((x2-x1)>0.0001){
132 ymed=ResolK1(xmed)-res;
137 ymed=Pol(xmed,2)-res;
140 if((y1*ymed)<0) x2=xmed;
141 if((y2*ymed)<0)x1=xmed;
142 if(ymed==0) return xmed;
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);
156 printf("k should be <=2\n");
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);
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);
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);
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();
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++;
197 //______________________________________________________________________
198 Int_t AliVertexingHFUtils::GetGeneratedMultiplicityInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
199 // counts generated particles in fgiven eta range
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++;
210 //______________________________________________________________________
211 Int_t AliVertexingHFUtils::GetGeneratedPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
212 // counts generated primary particles in given eta range
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++;
225 //______________________________________________________________________
226 Int_t AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray* arrayMC, Double_t mineta, Double_t maxeta){
227 // counts generated primary particles in given eta range
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++;
240 //______________________________________________________________________
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){
243 // Compute <pt> from 2D histogram M vs pt
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));
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);
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);
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);
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);
286 hMptBkgLo->Rebin(rebin);
287 hMptBkgHi->Rebin(rebin);
288 hMptSigReg->Rebin(rebin);
292 TH1F* hMptBkgLoScal=(TH1F*)hMptBkgLo->Clone("hPtBkgLoScalBin");
293 hMptBkgLoScal->Scale(bkgSig/bkgLow);
294 TH1F* hMptBkgHiScal=(TH1F*)hMptBkgHi->Clone("hPtBkgHiScalBin");
295 hMptBkgHiScal->Scale(bkgSig/bkgHi);
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.);
304 averagePt = hMptSig->GetMean();
305 errorPt = hMptSig->GetMeanError();
310 delete hMptBkgLoScal;
311 delete hMptBkgHiScal;
317 //____________________________________________________________________________
318 Double_t AliVertexingHFUtils::GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
319 // true impact parameter calculation for Dzero
321 if(!partD || !arrayMC || !mcHeader) return 99999.;
322 Int_t code=partD->GetPdgCode();
323 if(TMath::Abs(code)!=421) return 99999.;
326 mcHeader->GetVertex(vtxTrue);
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++){
337 Int_t nDau=partD->GetNDaughters();
338 Int_t labelFirstDau = partD->GetDaughter(0);
340 for(Int_t iDau=0; iDau<2; iDau++){
341 Int_t ind = labelFirstDau+iDau;
342 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
344 printf("Daughter particle not found in MC array");
347 pXdauTrue[iDau]=part->Px();
348 pYdauTrue[iDau]=part->Py();
349 pZdauTrue[iDau]=part->Pz();
355 Double_t d0dummy[2]={0.,0.};
356 AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
357 return aodDvsMC.ImpParXY();
360 //____________________________________________________________________________
361 Double_t AliVertexingHFUtils::GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
362 // true impact parameter calculation for Dplus
364 if(!partD || !arrayMC || !mcHeader) return 99999.;
365 Int_t code=partD->GetPdgCode();
366 if(TMath::Abs(code)!=411) return 99999.;
369 mcHeader->GetVertex(vtxTrue);
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++){
380 Int_t nDau=partD->GetNDaughters();
381 Int_t labelFirstDau = partD->GetDaughter(0);
383 for(Int_t iDau=0; iDau<3; iDau++){
384 Int_t ind = labelFirstDau+iDau;
385 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
387 printf("Daughter particle not found in MC array");
390 pXdauTrue[iDau]=part->Px();
391 pYdauTrue[iDau]=part->Py();
392 pZdauTrue[iDau]=part->Pz();
396 for(Int_t iDau=0; iDau<2; iDau++){
397 Int_t ind = labelFirstDau+iDau;
398 AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
400 printf("Daughter particle not found in MC array");
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();
410 Int_t nDauRes=part->GetNDaughters();
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));
417 printf("Daughter particle not found in MC array");
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();
433 printf("Wrong number of decay prongs");
438 Double_t d0dummy[3]={0.,0.,0.};
439 AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
440 return aodDvsMC.ImpParXY();
446 //____________________________________________________________________________
447 Double_t AliVertexingHFUtils::GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult) {
449 // Correct the number of accepted tracklets based on a TProfile Hist
453 if(TMath::Abs(vtxZ)>10.0){
454 // printf("ERROR: Z vertex out of range for correction of multiplicity\n");
455 return uncorrectedNacc;
459 printf("ERROR: Missing TProfile for correction of multiplicity\n");
460 return uncorrectedNacc;
463 Double_t localAvg = estimatorAvg->GetBinContent(estimatorAvg->FindBin(vtxZ));
465 Double_t deltaM = uncorrectedNacc*(refMult/localAvg - 1);
467 Double_t correctedNacc = uncorrectedNacc + (deltaM>0 ? 1 : -1) * gRandom->Poisson(TMath::Abs(deltaM));
469 return correctedNacc;
471 //______________________________________________________________________
472 TString AliVertexingHFUtils::GetGenerator(Int_t label, AliAODMCHeader* header){
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;