/************************************************************************* * Copyright(c) 1998-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. * **************************************************************************/ //////////////////////////////////////////////////////////////////////////// // Class to perform pt-spectra (and ptMin-spectra) extraction of mothers // // particles starting from measured pt-spectra of daughter particles // // that come from inclusive decays. // // E.g.: B->J/psi+X , B->e+X, B->D0+X, etc. // // // // In order to use this class, one first has to run a simulation // // (only kinematics) of the decay channel under study. The analysis // // can be runned using the class AliAnalysisTaskPtMothFromPtDaugh // // which loops on events to create a TNtupla that stores just // // kinematics informations for mothers and daughters of the decay // // under study (this is made in order to speed up). // // // // Therefore the standard inputs of this class are: // // (1) The TNtupla (created by the task using a TChain of galice.root) // // (2) pT-spectrum of the daughter particles // // // // Output would be the pT (and pTMin) spectrum of the mother, based // // on the correction factors computed from the Kinematics.root files // // // // // // Authors: Giuseppe E. Bruno & Fiorella Fionda // // (Giuseppe.Bruno@ba.infn.it) (Fiorella.Fionda@ba.infn.it) // //////////////////////////////////////////////////////////////////////////// #include "TH1F.h" #include "TNtuple.h" #include "TFile.h" #include "TParticle.h" #include "TArrayI.h" #include "AliPtMothFromPtDaugh.h" #include "AliStack.h" #include "AliLog.h" ClassImp(AliPtMothFromPtDaugh) //________________________________________________________________ AliPtMothFromPtDaugh::AliPtMothFromPtDaugh() : TNamed("AliPtMoth","AliPtMoth"), fDecayKine(0x0), fWij(0x0), fFi(0x0), fWijMin(0x0), fFiMin(0x0), fHistoPtDaughter(0x0), fHistoPtMothers(0x0), fHistoPtMinMothers(0x0), fMothers(0x0), fDaughter(0), fyMothMax(0), fyMothMin(0), fyDaughMax(0), fyDaughMin(0), fUseEta(kFALSE), fAnalysisMode(kUserAnalysis) { // // Default constructor // } //________________________________________________________________ AliPtMothFromPtDaugh::AliPtMothFromPtDaugh(const char* name, const char* title) : TNamed(name,title), fDecayKine(0x0), fWij(0x0), fFi(0x0), fWijMin(0x0), fFiMin(0x0), fHistoPtDaughter(0x0), fHistoPtMothers(0x0), fHistoPtMinMothers(0x0), fMothers(0x0), fDaughter(0), fyMothMax(0), fyMothMin(0), fyDaughMax(0), fyDaughMin(0), fUseEta(kFALSE), fAnalysisMode(kUserAnalysis) { // // Named constructor // } //________________________________________________________________ AliPtMothFromPtDaugh::~AliPtMothFromPtDaugh() { // // Default destructor // if(fDecayKine) {delete fDecayKine;} fDecayKine=0; if(fMothers) {delete fMothers;} fMothers=0; if(fHistoPtMothers) { delete fHistoPtMothers; } fHistoPtMothers=0; if(fHistoPtMinMothers) { delete fHistoPtMinMothers;} fHistoPtMinMothers=0; if(fHistoPtDaughter) { delete fHistoPtDaughter; fHistoPtDaughter=0;} } //______________________________________________________________________________________ AliPtMothFromPtDaugh::AliPtMothFromPtDaugh(const AliPtMothFromPtDaugh& extraction) : TNamed(extraction), fDecayKine(0), fWij(0), fFi(0), fWijMin(0), fFiMin(0), fHistoPtDaughter(0), fHistoPtMothers(0), fHistoPtMinMothers(0), fMothers(0), fDaughter(extraction.fDaughter), fyMothMax(extraction.fyMothMax), fyMothMin(extraction.fyMothMin), fyDaughMax(extraction.fyDaughMax), fyDaughMin(extraction.fyDaughMin), fUseEta(extraction.fUseEta), fAnalysisMode(extraction.fAnalysisMode) { // Copy constructor if(extraction.fHistoPtDaughter) fHistoPtDaughter =(TH1F*)extraction.fHistoPtDaughter->Clone("fHistoPtDaughter_copy"); if(extraction.fHistoPtMothers) fHistoPtMothers = (TH1F*)extraction.fHistoPtMothers->Clone("fHistoPtMothers_copy"); if(extraction.fHistoPtMinMothers) fHistoPtMinMothers = (TH1F*)extraction.fHistoPtMinMothers->Clone("fHistoPtMinMothers_copy"); if(extraction.fWij){ fWij = new Double_t*[2*(fHistoPtMothers->GetNbinsX())]; for(Int_t i=0;i<2*(fHistoPtMothers->GetNbinsX());i++) *(fWij+i)=new Double_t[fHistoPtDaughter->GetNbinsX()]; for(Int_t i=0;i<2*(fHistoPtMothers->GetNbinsX());i++){ for(Int_t j=0;jGetNbinsX();j++){ fWij[i][j]=extraction.fWij[i][j]; } } } if(extraction.fFi){ fFi=new Double_t[2*(fHistoPtMothers->GetNbinsX())]; for(Int_t i=0;i<2*(fHistoPtMothers->GetNbinsX());i++) fFi[i]=extraction.fFi[i]; } if(extraction.fWijMin){ fWijMin = new Double_t*[2*(fHistoPtMinMothers->GetNbinsX())]; for(Int_t i=0;i<2*(fHistoPtMinMothers->GetNbinsX());i++) *(fWijMin+i)=new Double_t[fHistoPtDaughter->GetNbinsX()]; for(Int_t i=0;i<2*(fHistoPtMinMothers->GetNbinsX());i++){ for(Int_t j=0;jGetNbinsX();j++){ fWijMin[i][j]=extraction.fWijMin[i][j]; } } } if(extraction.fFiMin){ fFiMin=new Double_t[2*(fHistoPtMinMothers->GetNbinsX())]; for(Int_t i=0;i<2*(fHistoPtMinMothers->GetNbinsX());i++) fFiMin[i]=extraction.fFiMin[i]; } if(extraction.fDecayKine) fDecayKine = (TNtuple*)extraction.fDecayKine->CloneTree(); if(extraction.fMothers) fMothers = new TArrayI(*(extraction.fMothers)); extraction.Copy(*this); } //______________________________________________________________________________________ AliPtMothFromPtDaugh& AliPtMothFromPtDaugh::operator=(const AliPtMothFromPtDaugh &extraction) { // operator assignment if (this!=&extraction) { this->~AliPtMothFromPtDaugh(); new(this) AliPtMothFromPtDaugh(extraction); } return *this; } //______________________________________________________________________________________ Bool_t AliPtMothFromPtDaugh::CreateWeights() { // Set mothers and daughters pdg codes if not // Put a control if mothers histograms binning-size, rapidity (or pseudoapidity) range // are setting. Read daughter histogram (histName) from the file (pathFileName) // Initialize dimensions for correction factors DeleteWeights(); if(!fMothers) {AliError("Set pdg codes of mothers by SetPdgMothers!\n"); return kFALSE;} if(!fDaughter) {AliError("Set pdg code of daughter by SetPdgDaughter!\n"); return kFALSE;} if(!fHistoPtDaughter) { AliError("Daughter histogram doesn't exist! \n"); return kFALSE;} //Set Rapidity or Pseudorapidity range for mothers if not if(!fyMothMax || !fyMothMin ){ AliError("Set rapidity range or pseudoRapidity range for mothers: use SetYmothers(ymin,ymax) or SetEtaMothers(etamin,etamax)"); return kFALSE;} if(!fyDaughMax || !fyDaughMin){ AliError("Set rapidity range or pseudoRapidity range for daughters:use SetYdaughters(ymin,ymax) or SetEtaDaughters(etamin,etamax)"); return kFALSE;} if(!fHistoPtMothers) {AliError("Call method SetBinsPtMoth to set pT-histogram "); return kFALSE;} if(!fHistoPtMinMothers){AliError("Call method SetBinsPtMinMoth to set pTmin-histogram "); return kFALSE;} Int_t nbinsM=(fHistoPtMothers->GetNbinsX()+2); Int_t nbinsD=(fHistoPtDaughter->GetNbinsX()+2); Int_t nbinsMmin=(fHistoPtMinMothers->GetNbinsX()+2); //Create pointers for weights to reconstruct daughter and mothers pT-spectra fWij=new Double_t*[2*nbinsM]; for(Int_t i=0;i<2*nbinsM;i++) {*(fWij+i)=new Double_t[nbinsD];} fFi=new Double_t[2*nbinsM]; fWijMin=new Double_t*[2*nbinsMmin]; for(Int_t i=0;i<2*nbinsMmin;i++) {*(fWijMin+i)=new Double_t[nbinsD];} fFiMin=new Double_t[2*nbinsMmin]; AliInfo(Form("Pt-mothers distribution: pt_min = %f pt_max=%f n_bins=%d \n", fHistoPtMothers->GetBinLowEdge(1),fHistoPtMothers->GetBinLowEdge(nbinsM-1),nbinsM-2)); AliInfo(Form("PtMinimum-mothers distribution: pt_min = %f pt_max=%f n_bins=%d \n", fHistoPtMinMothers->GetBinLowEdge(1),fHistoPtMinMothers->GetBinLowEdge(nbinsMmin-1),nbinsMmin-2)); AliInfo(Form("Pt-daughters distribution: pt_min = %f pt_max=%f n_bins=%d \n", fHistoPtDaughter->GetBinLowEdge(1),fHistoPtDaughter->GetBinLowEdge(nbinsD-1),nbinsD-2)); return kTRUE; } //______________________________________________________________________________________ void AliPtMothFromPtDaugh::DeleteWeights() { //delete correction factors //delete histogram of daughters if(!fHistoPtMothers || !fHistoPtMinMothers) {AliError("Mothers histograms don't exist! Cannot delete correction factors"); return;} Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; Int_t nbinsMmin=fHistoPtMinMothers->GetNbinsX()+2; if(fWij){ for(Int_t i=0; i<(2*nbinsM); i++) delete fWij[i]; delete [] fWij; fWij=0; } if(fFi) { delete fFi; fFi=0;} if(fWijMin){ for(Int_t i=0; i<(2*nbinsMmin); i++) delete fWijMin[i]; delete [] fWijMin; fWijMin=0; } if(fFiMin) { delete fFiMin; fFiMin=0;} return; } //______________________________________________________________________________________ Bool_t AliPtMothFromPtDaugh::ReadHistoPtDaught(const TH1F *hist) { //Initialize daughter histograms with hist if(!hist) {AliError("Set correct histograms of daughter! It doesn't exist!\n"); return kFALSE;} if(fHistoPtDaughter) delete fHistoPtDaughter; fHistoPtDaughter = (TH1F*)hist->Clone(); return kTRUE; } //______________________________________________________________________________________ Double_t* AliPtMothFromPtDaugh::SetBinsSize(Double_t ptmin, Double_t ptmax,Int_t nbins,Double_t alpha) { // return a pointer of double which contains the binning size for (mothers) histograms: // alpha have to be > 0 and: // alpha = 1 equal binning size // alpha < 1 increasing " // alpha > 1 decreasing " if(ptmin<0 || ptmax<0 || nbins<=0 || alpha<=0) {AliError("Set correct bin-size: ptmin>=0, ptmax>=0, nbins>0, alpha>0! \n"); return 0;} Double_t *edgebin = new Double_t[nbins+1]; Double_t ptmin1=TMath::Power(ptmin,alpha); Double_t ptmax1=TMath::Power(ptmax,alpha); Double_t size=(ptmax1-ptmin1)/nbins; for(Int_t i=0;i 0 and: // alpha = 1 equal binning size // alpha < 1 increasing " // alpha > 1 decreasing " Double_t* edges = SetBinsSize(ptmin,ptmax,nbins,alpha); SetBinsPtMoth(nbins,edges); delete edges; return; } //______________________________________________________________________________________ void AliPtMothFromPtDaugh::SetBinsPtMoth(Int_t nbins,Double_t *edgeBins) { //set bin size given by the pointer edgeBins for pt-spectrum of mothers: //the dimension of the pointer edgeBins is nbins+1 and the points //has to be written in increasing order if(nbins<0) {AliError("Numbers of bins should be > 0 !\n"); return;} if(fHistoPtMothers) {delete fHistoPtMothers; fHistoPtMothers=0;} fHistoPtMothers=new TH1F("fHistoPtMothers","Reconstructed p_{T}(Mothers)-spectrum",nbins,edgeBins); return; } //______________________________________________________________________________________ void AliPtMothFromPtDaugh::SetBinsPtMinMoth(Int_t nbins,Double_t *edgeBins) { //set bin size given by the pointer edgeBins for ptMin-spectrum of mothers: //the dimension of the pointer edgeBins is nbins+1 and the points //has to be written in increasing order if(nbins<0) {AliError("Numbers of bins should be > 0 !\n"); return;} if(fHistoPtMinMothers) {delete fHistoPtMinMothers; fHistoPtMinMothers=0;} fHistoPtMinMothers = new TH1F("fHistoPtMinMothers","Reconstructed p_{T}^{MIN}(Mothers)-spectrum",nbins,edgeBins); return; } //______________________________________________________________________________________ void AliPtMothFromPtDaugh::SetBinsPtMinMoth(Double_t ptmin, Double_t ptmax,Int_t nbins,Double_t alpha) { // Set bin size for ptMin-spectrum of mothers using SetBinsSize: // alpha have to be > 0 and: // alpha = 1 equal binning size // alpha < 1 increasing " // alpha > 1 decreasing " Double_t* edges = SetBinsSize(ptmin,ptmax,nbins,alpha); SetBinsPtMinMoth(nbins,edges); delete edges; return; } //______________________________________________________________________________________ void AliPtMothFromPtDaugh::SetPdgMothersPrivate(Int_t n_mothers,Int_t *pdgM) { // Set pdg codes of mothers given by the pointer pdgM for the analysis. // This is a private method. if(fMothers) { delete fMothers; fMothers = 0; } fMothers = new TArrayI(n_mothers,pdgM); return; } //______________________________________________________________________________________ void AliPtMothFromPtDaugh::SetPdgMothers(Int_t n_mothers,Int_t *pdgM) { // Set user pdg codes of mothers: first check // that the kUserAnalysis is the selected Analysis_mode. // If not print out a message of error. if(fAnalysisMode!=kUserAnalysis) { AliError("Nothing done: set the mode to kUserAnalysis first"); return; } //Set pdg codes of mothers given by the pointer pdgM for the analysis SetPdgMothersPrivate(n_mothers,pdgM); return; } //______________________________________________________________________________________ void AliPtMothFromPtDaugh::SetBeautyMothers() { // Set pdg codes of beauty particles: // B-mesons (1-24) B-barions (25-59) Int_t pdgBeauty[]={511,521,513,523,10511,10521,10513,10523,20513, 20523,515,525,531,10531,533,10533,205433,535,541,10541,543,10543, 20543,545,5122,5112,5212,5222,5114,5214,5224,5132,5232,5312,5322, 5314,5324,5332,5334,5142,5242,5412,5422,5414,5424,5342,5432,5434, 5442,5444,5512,5522,5514,5524,5532,5534,5542,5544,5554}; Int_t *pdgB=new Int_t[59]; for(Int_t i=0;i<59;i++) pdgB[i]=pdgBeauty[i]; SetPdgMothersPrivate(59,pdgB); } //______________________________________________________________________________________ void AliPtMothFromPtDaugh::InitDefaultAnalysis() { // Set mothers and daughter pdg codes depending from the selected analysis. // case kUserAnalysis: mothers and daughter are set by user (nothing to be done) // case kBtoJPSI: inclusive B-> J/Psi + X channels // case kBtoEle: inclusive B-> e + X channels // case kBtoMuon: inclusive B-> mu + X channels // case kBtoD0: inclusive B-> D0 + X channels switch(fAnalysisMode) { case kUserAnalysis: break; case kBtoJPSI: SetBeautyMothers(); fDaughter = 443; break; case kBtoEle: SetBeautyMothers(); fDaughter = 11; break; case kBtoMuon: SetBeautyMothers(); fDaughter = 13; break; case kBtoD0: SetBeautyMothers(); fDaughter = 421; break; } } //______________________________________________________________________________________ Double_t* AliPtMothFromPtDaugh::GetBinsSize(const TH1F *hist,Int_t &n) const { // return the binning size of the histogram hist // n return the number of bins Double_t* edges = (Double_t*)hist->GetXaxis()->GetXbins()->GetArray(); n=hist->GetNbinsX(); return edges; } //______________________________________________________________________________________ Bool_t AliPtMothFromPtDaugh::GetEtaMothers(Double_t &etaMin, Double_t &etaMax) const { // method to get the bounds of the pseudorapidity range // for mothers. Return kTRUE if pseudorapidity is used and put // pseudorapidity edges in etaMin and etaMax if(fUseEta == kFALSE){ AliError("You are using RAPIDITY range! \n"); etaMin = 0.; etaMax = 0.; return kFALSE; } etaMin = fyMothMin; etaMax = fyMothMax; return kTRUE; } //______________________________________________________________________________________ Bool_t AliPtMothFromPtDaugh::GetEtaDaughter(Double_t &etaMin, Double_t &etaMax) const { // method to get the bounds of the pseudorapidity range // for daughters. Return kTRUE if pseudorapidity is used and put // pseudorapidity edges in etaMin and etaMax if(fUseEta == kFALSE){ AliError("You are using RAPIDITY range! \n"); etaMin = 0.; etaMax = 0.; return kFALSE; } etaMin = fyDaughMin; etaMax = fyDaughMax; return kTRUE; } //______________________________________________________________________________________ Bool_t AliPtMothFromPtDaugh::GetYMothers(Double_t &yMin, Double_t &yMax) const { // method to get the bounds of the rapidity range // for mothers. Return kTRUE if rapidity is used and put // rapidity edges in yMin and yMax if(fUseEta == kTRUE){ AliError("You are using PSEUDORAPIDITY range! \n"); yMin = 0.; yMax = 0.; return kFALSE; } yMin = fyMothMin; yMax = fyMothMax; return kTRUE; } //______________________________________________________________________________________ Bool_t AliPtMothFromPtDaugh::GetYDaughter(Double_t &yMin, Double_t &yMax) const { // method to get the bounds of the rapidity range // for daughters. Return kTRUE if rapidity is used and put // rapidity edges in yMin and yMax if(fUseEta == kTRUE){ AliError("You are using PSEUDORAPIDITY range! \n"); yMin = 0.; yMax = 0.; return kFALSE; } yMin = fyDaughMin; yMax = fyDaughMax; return kTRUE; } //______________________________________________________________________________________ void AliPtMothFromPtDaugh::SetPdgDaugh(Int_t pdgD) { // Set pdg code for daughter particle. Check // that the kUserAnalysis is the selected Analysis_mode. // If not print out a message of error. switch(fAnalysisMode) { case kUserAnalysis: fDaughter = pdgD; break; case kBtoJPSI: if(pdgD!= 443) {AliError("Nothing done: first change AnalysisMode to kUserAnalysis");} else {fDaughter = pdgD;} break; case kBtoEle: if(pdgD!= 11) {AliError("Nothing done: first change AnalysisMode to kUserAnalysis");} else {fDaughter = pdgD;} break; case kBtoMuon: if(pdgD!= 13) {AliError("Nothing done: first change AnalysisMode to kUserAnalysis");} else {fDaughter = pdgD;} break; case kBtoD0: if(pdgD!= 421) {AliError("Nothing done: first change AnalysisMode to kUserAnalysis");} else {fDaughter = pdgD;} break; } return; } //______________________________________________________________________________________ Int_t* AliPtMothFromPtDaugh::GetPdgMothers(Int_t &n_mothers) const { // return the pointer to the array of pdg codes of mothers particles // if it exist. Put its dimension in n_mothers if(!fMothers) {AliError("Mothers pdg are not defined! \n"); return 0x0;} n_mothers = fMothers->GetSize(); return fMothers->GetArray(); } //______________________________________________________________________________________ Double_t AliPtMothFromPtDaugh::GetW(Int_t i,Int_t j) const { // Return value of correction factors Wij at the position i (pt-mothers bin index)- // j (pt-daughter bin index). Bin 0 is the underflow, bin nbins+1 the overflow. // If Wij don't exist or the indices i or j are out of the variability range return 0 if(!fWij) {AliError("Correction factors Wij are not been evaluated yet!\n"); return 0;} if(!fHistoPtMothers) {AliError("mothers pt-histogram doesn't exist!\n"); return 0;} if(!fHistoPtDaughter) {AliError("daughters pt-histogram doesn't exist!\n"); return 0;} Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; Int_t nbinsD=fHistoPtDaughter->GetNbinsX()+2; if(i<0 || i>nbinsM-1) {AliError(Form("Index i out of range: %d =< i =< %d. Bin 0 is the underflow - bin %d is the overflow ",0,nbinsM-1,nbinsM-1)); return 0;} if(j<0 || j>nbinsD-1) {AliError(Form("Index j out of range: %d =< j =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsD-1,nbinsD-1)); return 0;} return fWij[i][j]; } //______________________________________________________________________________________ Double_t AliPtMothFromPtDaugh::GetStatErrW(Int_t i,Int_t j) const { // Return value of statistical error on correction factors Wij at the position // i (pt-mothers bin index)- j (pt-daughter bin index). Bin 0 is the underflow, // bin nbins+1 the overflow. If Wij don't exist or the indices i or j are out of the // variability range return 0 if(!fHistoPtMothers) {AliError("mothers pt-histogram doesn't exist!\n"); return 0;} if(!fHistoPtDaughter) {AliError("daughters pt-histogram doesn't exist!\n"); return 0;} if(!fWij) {AliError("Correction factors Wij are not been evaluated yet!\n"); return 0;} Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; Int_t nbinsD=fHistoPtDaughter->GetNbinsX()+2; if(i<0 || i>nbinsM-1) {AliError(Form("Index i out of range: %d =< i =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsM-1,nbinsM-1)); return 0;} if(j<0 || j>nbinsD-1) {AliError(Form("Index j out of range: %d =< j =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsD-1,nbinsD-1)); return 0;} return fWij[i+nbinsM][j]; } //______________________________________________________________________________________ Double_t AliPtMothFromPtDaugh::GetF(Int_t i) const { // Return value of correction factors Fi at the position i (pt-mothers bin index). // Bin 0 is the underflow, bin nbins+1 the overflow. // If Fi don't exist or the index i is out of the variability range return 0 if(!fFi) {AliError("Correction factors Fi are not been evaluated yet!\n"); return 0;} if(!fHistoPtMothers) {AliError("mothers pt-histogram doesn't exist!\n"); return 0;} Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; if(i<0 || i>nbinsM-1) {AliError(Form("Index i out of range: %d =< i =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsM-1,nbinsM-1)); return 0;} return fFi[i]; } //______________________________________________________________________________________ Double_t AliPtMothFromPtDaugh::GetStatErrF(Int_t i) const { // Return statistical error on correction factors Fi at the position i (pt-mothers bin index). // Bin 0 is the underflow, bin nbins+1 the overflow. // If Fi don't exist or the index i is out of the variability range return 0 if(!fFi) {AliError("Correction factors Fi are not been evaluated yet!\n"); return 0;} if(!fHistoPtMothers) {AliError("mothers pt-histogram doesn't exist!\n"); return 0;} Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; if(i<0 || i>nbinsM-1) {AliError(Form("Index i out of range: %d =< i =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsM-1,nbinsM-1)); return 0;} return fFi[i+nbinsM]; } //______________________________________________________________________________________ Double_t AliPtMothFromPtDaugh::GetWmin(Int_t i,Int_t j) const { // Return value of correction factors Wij_min at the position i (ptMin-mothers bin index)- // j (pt-daughter bin index). Bin 0 is the underflow, bin nbins+1 the overflow. // If Wij_min don't exist or the indices i or j are out of the variability range return 0 if(!fWijMin) {AliError("Correction factors Wij_min are not been evaluated yet!\n"); return 0;} if(!fHistoPtMinMothers) {AliError("mothers ptMin-histogram doesn't exist!\n"); return 0;} if(!fHistoPtDaughter) {AliError("daughters pt-histogram doesn't exist!\n"); return 0;} Int_t nbinsMmin=fHistoPtMinMothers->GetNbinsX()+2; Int_t nbinsD=fHistoPtDaughter->GetNbinsX()+2; if(i<0 || i>nbinsMmin-1) {AliError(Form("Index i out of range: %d =< i =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsMmin-1,nbinsMmin-1)); return 0;} if(j<0 || j>nbinsD-1) {AliError(Form("Index j out of range: %d =< j =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsD-1,nbinsD-1)); return 0;} return fWijMin[i][j]; } //______________________________________________________________________________________ Double_t AliPtMothFromPtDaugh::GetStatErrWmin(Int_t i,Int_t j) const { // Return value of statistical error on correction factors Wij_min at the position // i (ptMin-mothers bin index)- j (pt-daughter bin index). Bin 0 is the underflow, // bin nbins+1 the overflow. If Wij_min don't exist or the indices i or j are out of the // variability range return 0 if(!fWijMin) {AliError("Correction factors Wij_min are not been evaluated yet!\n"); return 0;} if(!fHistoPtMinMothers) {AliError("mothers ptMin-histogram doesn't exist!\n"); return 0;} if(!fHistoPtDaughter) {AliError("daughters pt-histogram doesn't exist!\n"); return 0;} Int_t nbinsMmin=fHistoPtMinMothers->GetNbinsX()+2; Int_t nbinsD=fHistoPtDaughter->GetNbinsX()+2; if(i<0 || i>nbinsMmin-1) {AliError(Form("Index i out of range: %d =< i =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsMmin-1,nbinsMmin-1)); return 0;} if(j<0 || j>nbinsD-1) {AliError(Form("Index j out of range: %d =< j =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsD-1,nbinsD-1)); return 0;} return fWijMin[i+nbinsMmin][j]; } //______________________________________________________________________________________ Double_t AliPtMothFromPtDaugh::GetFmin(Int_t i) const { // Return value of correction factors Fi_min at the position i (ptMin-mothers bin index). // Bin 0 is the underflow, bin nbins+1 the overflow. // If Fi_min don't exist or the index i is out of the variability range return 0 if(!fFiMin) {AliError("Correction factors Fi_min are not been evaluated yet!\n"); return 0;} if(!fHistoPtMinMothers) {AliError("mothers ptMin-histogram doesn't exist!\n"); return 0;} Int_t nbinsMmin=fHistoPtMothers->GetNbinsX()+2; if(i<0 || i>nbinsMmin-1) {AliError(Form("Index i out of range: %d =< i =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsMmin-1,nbinsMmin-1)); return 0;} return fFiMin[i]; } //______________________________________________________________________________________ Double_t AliPtMothFromPtDaugh::GetStatErrFmin(Int_t i) const { // Return statistical error on correction factors Fi_min at the position i (ptMin-mothers bin index). // Bin 0 is the underflow, bin nbins+1 the overflow. // If Fi_min don't exist or the index i is out of the variability range return 0 if(!fFiMin) {AliError("Correction factors Fi_min are not been evaluated yet!\n"); return 0;} if(!fHistoPtMinMothers) {AliError("mothers ptMin-histogram doesn't exist!\n"); return 0;} Int_t nbinsMmin=fHistoPtMothers->GetNbinsX()+2; if(i<0 || i>nbinsMmin-1) {AliError(Form("Index i out of range: %d =< i =< %d. Bin 0 is the underflow - bin %d is the overflow.",0,nbinsMmin-1,nbinsMmin-1)); return 0;} return fFiMin[i+nbinsMmin]; } //______________________________________________________________________________________ Bool_t AliPtMothFromPtDaugh::IsMothers(Int_t pdgCode) { //return kTRUE if pdgCode is in the list of pdg codes of mothers Int_t dim = fMothers->GetSize(); for(Int_t i=0;iGetAt(i); if(pdgCode == pdgMother) return kTRUE; } return kFALSE; } //______________________________________________________________________________________ Bool_t AliPtMothFromPtDaugh::IsSelectedDaugh(const TParticle *part, Int_t &labelDaugh, AliStack *const stack) { // return kTRUE if particle part has the selected daughter // if yes put the label of the track in labelDaugh TParticle *daugh; Int_t nDg=part->GetNDaughters(); if(nDg<=0) {AliError("I have no daugh!\n"); return kFALSE;} for(Int_t i=part->GetFirstDaughter();iGetLastDaughter()+1;i++) { daugh=stack->Particle(i); if(TMath::Abs(daugh->GetPdgCode())==fDaughter) { labelDaugh=i; return kTRUE; } } return kFALSE; } //______________________________________________________________________________________ Bool_t AliPtMothFromPtDaugh::Rapidity(const TParticle *particle, Double_t &y) { // Evaluated rapidity of particle and put it in y. Return kFALSE if // cannot compute rapidity y=-999; if(particle->Energy()-particle->Pz()<=0) return kFALSE; y=0.5*TMath::Log((particle->Energy()+particle->Pz())/(particle->Energy()-particle->Pz())); return kTRUE; } //______________________________________________________________________________________ Int_t AliPtMothFromPtDaugh::GiveBinIndex(Double_t Ptpart,const TH1F *ptHist) const { // Return the bin index of pt respect to the binning size of ptHist // bin 0 is the underflow - nbins+1 is the overflow Int_t nbins=ptHist->GetNbinsX(); Int_t index=0; for(Int_t i=1;iGetBinLowEdge(i))) { index=i-1; break; } } return index; } //______________________________________________________________________________________ Bool_t AliPtMothFromPtDaugh::CutDaugh(Double_t yD,Double_t ptD) { // put a control for rapidity yD and transverse momentum ptD of daughter // return kTRUE if fyDaughMin < yD < fyDaughMax and ptMinDaugh < ptD < ptMaxDaugh Double_t ptMinDaugh = fHistoPtDaughter->GetBinLowEdge(1); Double_t ptMaxDaugh = fHistoPtDaughter->GetBinLowEdge(fHistoPtDaughter->GetNbinsX()); if( yD < fyDaughMin || yD > fyDaughMax ) return kFALSE; if( ptD > ptMaxDaugh || ptD < ptMinDaugh ) return kFALSE; return kTRUE; } //______________________________________________________________________________________ Bool_t AliPtMothFromPtDaugh::EvaluateWij() { // Evaluate correction factors using to extract the ptRaw and // ptMinRaw distributions. Statistical errors on those are computed too if(!fHistoPtMothers || !fHistoPtDaughter || !fHistoPtMinMothers) {AliError("Control mother and daughter histograms!\n"); return kFALSE;} Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; Int_t nbinsD=fHistoPtDaughter->GetNbinsX()+2; Int_t nbinsMmin=fHistoPtMinMothers->GetNbinsX()+2; Float_t pdgM, pdgD, pxM,pyM,pzM,yM,pxD,pyD,pzD,yD,etaM,etaD,cutVarM,cutVarD; Double_t* entries=new Double_t[nbinsD]; for(Int_t ii=0;iiSetBranchAddress("pdgM",&pdgM); fDecayKine->SetBranchAddress("pxM",&pxM); fDecayKine->SetBranchAddress("pyM",&pyM); fDecayKine->SetBranchAddress("pzM",&pzM); fDecayKine->SetBranchAddress("yM",&yM); fDecayKine->SetBranchAddress("etaM",&etaM); fDecayKine->SetBranchAddress("pdgD",&pdgD); fDecayKine->SetBranchAddress("pxD",&pxD); fDecayKine->SetBranchAddress("pyD",&pyD); fDecayKine->SetBranchAddress("pzD",&pzD); fDecayKine->SetBranchAddress("yD",&yD); fDecayKine->SetBranchAddress("etaD",&etaD); Double_t ptD,ptM; // Initialize correction factors for pT and pTmin if those exist if(!fWij) {AliError("Correction factors Wij have not been created!\n"); return kFALSE;} if(!fWijMin) {AliError("Correction factors Wij_min have not been created!\n"); return kFALSE;} for(Int_t ii=0;ii<(2*nbinsM);ii++){ for(Int_t jj=0;jjGetEntries(); Int_t fNcurrent=0; Int_t nb = (Int_t)fDecayKine->GetEvent(fNcurrent); for (Int_t iev=0; ievGetEvent(fNcurrent); continue;} if(cutVarM>fyMothMin && cutVarMGetEvent(fNcurrent); } for(Int_t jj=0;jj0){ fWij[ii][jj]=fWij[ii][jj]/entries[jj]; // evaluate statistical errors on fWij fWij[ii+nbinsM][jj]=(TMath::Sqrt(fWij[ii][jj]*(1-(fWij[ii][jj]/entries[jj]))))/entries[jj]; } else{ // if there are no entries in the bin-j of daughter distribution // set factor = -1 and error = 999 fWij[ii][jj]=-1; fWij[ii+nbinsM][jj]=999; } } for(Int_t ii=0;ii0){ fWijMin[ii][jj]=fWijMin[ii][jj]/entries[jj]; //evaluate statistical errors on fWijMin fWijMin[ii+nbinsMmin][jj] = (TMath::Sqrt(fWijMin[ii][jj]*(1-(fWijMin[ii][jj]/entries[jj]))))/entries[jj]; } else{ //if there are no entries set correction factor = -1 and error = -999 fWijMin[ii][jj]=-1; fWijMin[ii+nbinsMmin][jj]=999; } } } delete entries; return kTRUE; } //______________________________________________________________________________________ Bool_t AliPtMothFromPtDaugh::EvaluateFi() { // Evaluate acceptance correction factors that are applied on the // raw distributions. Statistical errors on those are computed too if(!fHistoPtMothers || !fHistoPtMinMothers) {AliError("Control mother histograms!\n"); return kFALSE;} Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; Int_t nbinsMmin=fHistoPtMinMothers->GetNbinsX()+2; Double_t* entries=new Double_t[nbinsM]; Double_t* entries1=new Double_t[nbinsMmin]; if(!fFi) {AliError("Correction factors Fi have not been created!\n"); return kFALSE;} if(!fFiMin) {AliError("Correction factors Fi_min have not been created!\n"); return kFALSE;} //initialize the correction factor for pT and pTmin for(Int_t ii=0;iiSetBranchAddress("pdgM",&pdgM); fDecayKine->SetBranchAddress("pxM",&pxM); fDecayKine->SetBranchAddress("pyM",&pyM); fDecayKine->SetBranchAddress("pzM",&pzM); fDecayKine->SetBranchAddress("yM",&yM); fDecayKine->SetBranchAddress("etaM",&etaM); fDecayKine->SetBranchAddress("pdgD",&pdgD); fDecayKine->SetBranchAddress("pxD",&pxD); fDecayKine->SetBranchAddress("pyD",&pyD); fDecayKine->SetBranchAddress("pzD",&pzD); fDecayKine->SetBranchAddress("yD",&yD); fDecayKine->SetBranchAddress("etaD",&etaD); Double_t ptD,ptM; Int_t nentries = (Int_t)fDecayKine->GetEntries(); Int_t fNcurrent=0; Int_t nb = (Int_t)fDecayKine->GetEvent(fNcurrent); for (Int_t iev=0; ievfyMothMax){ fNcurrent++; nb = (Int_t)fDecayKine->GetEvent(fNcurrent); continue;} entries[i]++; for(Int_t k=0; kGetEvent(fNcurrent); continue;} fFi[i]+=1.; for(Int_t k=0; kGetEvent(fNcurrent); } for(Int_t ii=0;ii0){ fFi[ii]/=entries[ii]; fFi[ii+nbinsM]=(TMath::Sqrt(fFi[ii]*(1-(fFi[ii]/entries[ii]))))/entries[ii]; } else{ fFi[ii]=-1; fFi[ii+nbinsM]=999; } } for(Int_t ii=0;ii0){ fFiMin[ii]/=entries1[ii]; fFiMin[ii+nbinsMmin]=(TMath::Sqrt(fFiMin[ii]*(1-(fFiMin[ii]/entries1[ii]))))/entries1[ii]; } else {fFiMin[ii]=-1; fFiMin[ii+nbinsMmin]=999;} } delete entries; delete entries1; return kTRUE; } //______________________________________________________________________________________ Bool_t AliPtMothFromPtDaugh::EvaluatePtMothRaw(TH1F *histoPt, TH1F *histoPtMin) { //Apply the fWij and fWijMin on the daughter distribution //in order to evaluate the pt and ptMin raw distributions for mothers if(!fHistoPtMothers || !fHistoPtDaughter || !fHistoPtMinMothers) {AliError("Control mother and daughter histograms!\n"); return kFALSE;} Int_t nbinsD=fHistoPtDaughter->GetNbinsX()+2; Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; Int_t nbinsMmin=fHistoPtMinMothers->GetNbinsX()+2; Double_t lowedge=0.; Double_t wfill=0.; Double_t wMinfill=0.; for(Int_t j=0;jGetBinCenter(i); if(fWij[i][j]>=0) wfill=fWij[i][j]*(fHistoPtDaughter->GetBinContent(j)); histoPt->Fill(lowedge,wfill); } for(Int_t i=0;iGetBinLowEdge(i); if(fWijMin[i][j]>=0) wMinfill=fWijMin[i][j]*(fHistoPtDaughter->GetBinContent(j)); histoPtMin->Fill(lowedge,wMinfill); } } return kTRUE; } //______________________________________________________________________________________ Bool_t AliPtMothFromPtDaugh::EvaluateErrPt(Double_t *erStat) { // Evaluate the statistical error on the pt-mothers distribution. // sigmaX: contribution that came from the measured distibution // sigmaWij: contribution that came from the fWij factors // sigmaFi: contribution that came from the fFi factors if(!fHistoPtMothers || !fHistoPtDaughter) {AliError("Control mother(pt) and daughter histograms!\n"); return kFALSE;} Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; Int_t nbinsD=fHistoPtDaughter->GetNbinsX()+2; Double_t m=0; Double_t sigmaX, sigmaWij, sigmaFi; for(Int_t i=0;i=0){ sigmaX+=(((fWij[i][j])*(fWij[i][j]))*((fHistoPtDaughter->GetBinError(j))*(fHistoPtDaughter->GetBinError(j)))); sigmaWij+=((fHistoPtDaughter->GetBinContent(j))*(fHistoPtDaughter->GetBinContent(j)))*(fWij[i+nbinsM][j]*fWij[i+nbinsM][j]); sigmaFi+=(fWij[i][j])*(fHistoPtDaughter->GetBinContent(j)); } } if(fFi[i]>0) sigmaFi=((sigmaFi*sigmaFi)*(fFi[i+nbinsM]*fFi[i+nbinsM]))/(fFi[i]*fFi[i]); m=TMath::Sqrt(sigmaX+sigmaWij+sigmaFi); if(fFi[i]>0) erStat[i]=(1/fFi[i])*m; else erStat[i]=999; } return kTRUE; } //______________________________________________________________________________________ Bool_t AliPtMothFromPtDaugh::EvaluateErrPtMin(Double_t *erStat) { // Evaluate statistical error on ptMin mothers distribution // sigmaMinX: contribution that came from the measured distibution // sigmaMinWij: contribution that came from the fWijMin factors // sigmaMinFi: contribution that came from the fFiMin factors if(!fHistoPtDaughter || !fHistoPtMinMothers) {AliError("Control mother(ptMin) and daughter histograms!\n"); return kFALSE;} Int_t nbinsD=fHistoPtDaughter->GetNbinsX()+2; Int_t nbinsMmin=fHistoPtMinMothers->GetNbinsX()+2; Double_t m1=0; Double_t sigmaMinX; Double_t sigmaMinWij; Double_t sigmaMinFi; for(Int_t i=0;i=0){ sigmaMinX+=(((fWijMin[i][j])*(fWijMin[i][j]))*((fHistoPtDaughter->GetBinError(j))*(fHistoPtDaughter->GetBinError(j)))); sigmaMinWij+=((fHistoPtDaughter->GetBinContent(j))*(fHistoPtDaughter->GetBinContent(j)))*(fWijMin[i+nbinsMmin][j]*fWijMin[i+nbinsMmin][j]); sigmaMinFi+=(fWijMin[i][j])*(fHistoPtDaughter->GetBinContent(j)); } } if(fFiMin[i]>0) sigmaMinFi=((sigmaMinFi*sigmaMinFi)*(fFiMin[i+nbinsMmin]*fFiMin[i+nbinsMmin]))/(fFiMin[i]*fFiMin[i]); m1=TMath::Sqrt(sigmaMinX+sigmaMinWij+sigmaMinFi); if(fFiMin[i]>0) erStat[i]=(1/fFiMin[i])*m1; else erStat[i]=999; } return kTRUE; } //______________________________________________________________________________________ Bool_t AliPtMothFromPtDaugh::EvaluatePtMoth() { // Evaluate pt and ptMin distribution for mothers // First evaluate the sigma raw distribution by calling EvaluatePtMothRaw // then evaluate the pt and ptMin mothers distribution. // Statistical errors on those distributions are evaluated too. if(!EvaluateWij()) return kFALSE; if(!EvaluateFi()) return kFALSE; // reset pt and ptMin mothers histograms fHistoPtMothers->Reset(); fHistoPtMinMothers->Reset(); TH1F *histoPt=(TH1F*)fHistoPtMothers->Clone(); TH1F *histoPtMin=(TH1F*)fHistoPtMinMothers->Clone(); EvaluatePtMothRaw(histoPt,histoPtMin); Int_t nbinsM=fHistoPtMothers->GetNbinsX()+2; Int_t nbinsMmin=fHistoPtMinMothers->GetNbinsX()+2; Double_t *erPtStat=new Double_t[nbinsM+2]; EvaluateErrPt(erPtStat); Double_t *erPtMinStat=new Double_t[nbinsMmin+2]; EvaluateErrPtMin(erPtMinStat); Double_t lowedge=0; Double_t fwfill; Double_t fMinfill; for(Int_t i=0;iGetBinCenter(i); if(fFi[i]>0){ fwfill=(histoPt->GetBinContent(i))/fFi[i]; fHistoPtMothers->Fill(lowedge,fwfill); fHistoPtMothers->SetBinError(i,erPtStat[i]); } } for(Int_t i=0;iGetBinCenter(i); if(fFiMin[i]>0){ fMinfill=(histoPtMin->GetBinContent(i))/fFiMin[i]; fHistoPtMinMothers->Fill(lowedge,fMinfill); fHistoPtMinMothers->SetBinError(i,erPtMinStat[i]); } } return kTRUE; } //______________________________________________________________________________________ void AliPtMothFromPtDaugh::WritePtMothHistoToFile(char *fileOutName) { // Write pt and ptMin histograms of mothers in a file // with name fileOutName. Default name is "Mothers.root". AliError(Form("Write mothers histograms in the file %s \n",fileOutName)); if(!fHistoPtMothers) {AliError("Cannot write pt-mothers histogram! It doesn't exists!"); return;} if(!fHistoPtMinMothers) { AliError("Cannot write ptMin-mothers histogram! It doesn't exists!"); return;} TFile *outFile = TFile::Open(fileOutName,"RECREATE"); outFile->cd(); fHistoPtMothers->Write(); fHistoPtMinMothers->Write(); outFile->Close(); return; }