[u/mrichter/AliRoot.git] / ANALYSIS / AliAnalysisTaskPIDResponse.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, 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.                  *
 **************************************************************************/

/* $Id: AliAnalysisTaskPIDResponse.cxx 43811 2010-09-23 14:13:31Z wiechula $ */
#include <TList.h>
#include <TVectorD.h>
#include <TObjArray.h>
#include <TH2.h>
#include <TFile.h>
#include <TPRegexp.h>
#include <TChain.h>

#include <AliAnalysisManager.h>
#include <AliInputEventHandler.h>
#include <AliVEventHandler.h>
#include <AliVEvent.h>
#include <AliVParticle.h>
#include <AliVTrack.h>
#include <AliLog.h>
#include <AliPID.h>
#include <AliPIDResponse.h>
#include <AliITSPIDResponse.h>
#include <AliTPCPIDResponse.h>
#include <AliTRDPIDResponse.h>
#include <AliTOFPIDResponse.h>

#include "AliAnalysisTaskPIDResponse.h"


ClassImp(AliAnalysisTaskPIDResponse)

//______________________________________________________________________________
AliAnalysisTaskPIDResponse::AliAnalysisTaskPIDResponse():
AliAnalysisTaskSE(),
fIsMC(kFALSE),
fTOFTimeZeroTypeUser(-1),
fTOFTimeZeroType(AliPIDResponse::kBest_T0),
fTOFres(100.),
fPIDResponse(0x0),
fListQA(0x0),
fListQAits(0x0),
fListQAtpc(0x0),
fListQAtrd(0x0),
fListQAtof(0x0),
fBeamType("PP"),
fLHCperiod(),
fMCperiodTPC(),
fRecoPass(0),
fRun(0),
fOldRun(0),
fArrPidResponseMaster(0x0)
{
  //
  // Dummy constructor
  //
}

//______________________________________________________________________________
AliAnalysisTaskPIDResponse::AliAnalysisTaskPIDResponse(const char* name):
AliAnalysisTaskSE(name),
fIsMC(kFALSE),
fTOFTimeZeroTypeUser(-1),
fTOFTimeZeroType(AliPIDResponse::kBest_T0),
fTOFres(100.),
fPIDResponse(0x0),
fListQA(0x0),
fListQAits(0x0),
fListQAtpc(0x0),
fListQAtrd(0x0),
fListQAtof(0x0),
fBeamType("PP"),
fLHCperiod(),
fMCperiodTPC(),
fRecoPass(0),
fRun(0),
fOldRun(0),
fArrPidResponseMaster(0x0)
{
  //
  // Default constructor
  //
  DefineInput(0,TChain::Class());
  DefineOutput(1,TList::Class());
}

//______________________________________________________________________________
AliAnalysisTaskPIDResponse::~AliAnalysisTaskPIDResponse()
{
  //
  // Destructor
  //

  delete fArrPidResponseMaster;
}

//______________________________________________________________________________
void AliAnalysisTaskPIDResponse::ExecNewRun()
{
  //
  // Things to Execute upon a new run 
  //

  SetRecoInfo();
    
  SetITSParametrisation();

  SetTPCPidResponseMaster();
  SetTPCParametrisation();

  fPIDResponse->GetTOFResponse().SetTimeResolution(fTOFres);
}

//______________________________________________________________________________
void AliAnalysisTaskPIDResponse::UserCreateOutputObjects()
{
  //
  // Create the output QA objects
  //

  AliLog::SetClassDebugLevel("AliAnalysisTaskPIDResponse",10);

  //input hander
  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
  AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
  if (!inputHandler) AliFatal("Input handler needed");

  //pid response object
  inputHandler->CreatePIDResponse(fIsMC);
  fPIDResponse=inputHandler->GetPIDResponse();
  
  //
  fListQA=new TList;
  fListQAits=new TList;
  fListQAits->SetName("ITS");
  fListQAtpc=new TList;
  fListQAtpc->SetName("TPC");
  fListQAtrd=new TList;
  fListQAtrd->SetName("TRD");
  fListQAtof=new TList;
  fListQAtof->SetName("TOF");
  
  fListQA->Add(fListQAits);
  fListQA->Add(fListQAtpc);
  fListQA->Add(fListQAtrd);
  fListQA->Add(fListQAtof);

  SetupTTSqa();
  SetupTPCqa();
  SetupTRDqa();
  SetupTOFqa();

  PostData(1,fListQA);
}


//______________________________________________________________________________
void AliAnalysisTaskPIDResponse::UserExec(Option_t */*option*/)
{
  //
  // Setup the PID response functions and fill the QA histograms
  //

  AliVEvent *event=InputEvent();
  if (!event) return;
  fRun=event->GetRunNumber();

  if (fRun!=fOldRun){
    ExecNewRun();
    fOldRun=fRun;
  }

  Double_t timeZeroType=fTOFTimeZeroTypeUser;
  if (timeZeroType<0) timeZeroType=fTOFTimeZeroType;
  fPIDResponse->SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)timeZeroType);
  
  FillITSqa();
  FillTPCqa();
  FillTOFqa();

  PostData(1,fListQA);
}


//______________________________________________________________________________
void AliAnalysisTaskPIDResponse::SetRecoInfo()
{
  //
  // Set reconstruction information
  //
  
  //reset information
  fRecoPass=0;
  fLHCperiod="";
  fMCperiodTPC="";

  fBeamType="";
  
  //Get the current file to check the reconstruction pass (UGLY, but not stored in ESD... )
  AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
  AliVEventHandler *inputHandler=mgr->GetInputEventHandler();
  if (!inputHandler) return;
  
  TTree *tree= (TTree*)inputHandler->GetTree();
  TFile *file= (TFile*)tree->GetCurrentFile();
  
  if (!file) {
    AliError("Current file not found, cannot set reconstruction information");
    return;
  }

  fBeamType="PP";
  
  //find the period by run number (UGLY, but not stored in ESD and AOD... )
  if (fRun>=114737&&fRun<=117223)      { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1";  }
  else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1";  }
  else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
  else if (fRun>=127719&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
  else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
  else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
  else if (fRun>=136851&&fRun<=139517) { fLHCperiod="LHC10H"; fMCperiodTPC="LHC10H8"; fBeamType="PBPB"; }
  else if (fRun>=139699) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
  
  //find pass from file name (UGLY, but not stored in ESD... )
  TString fileName(file->GetName());
  if (fileName.Contains("/pass1")) {
    fRecoPass=1;
  } else if (fileName.Contains("/pass2")) {
    fRecoPass=2;
  } else if (fileName.Contains("/pass3")) {
    fRecoPass=3;
  }
  
}

//______________________________________________________________________________
void AliAnalysisTaskPIDResponse::SetITSParametrisation()
{
  //
  // Set the ITS parametrisation
  //
}

//______________________________________________________________________________
void AliAnalysisTaskPIDResponse::SetTPCPidResponseMaster()
{
  //
  // Load the TPC pid response functions from the OADB
  //

  //reset the PID response functions
  delete fArrPidResponseMaster;
  fArrPidResponseMaster=0x0;
  
  TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", AliAnalysisManager::GetOADBPath()));

  TFile f(fileName.Data());
  if (f.IsOpen() && !f.IsZombie()){
    fArrPidResponseMaster=dynamic_cast<TObjArray*>(f.Get("TPCPIDResponse"));
    fArrPidResponseMaster->SetOwner();
    f.Close();
  }

  if (!fArrPidResponseMaster){
    AliFatal("Could not retrieve the TPC pid response");
    return;
  }
}

//______________________________________________________________________________
void AliAnalysisTaskPIDResponse::SetTPCParametrisation()
{
  //
  // Change BB parametrisation for current run
  //
  
  if (fLHCperiod.IsNull()) {
    AliFatal("No period set, not changing parametrisation");
    return;
  }
  
  //
  // Set default parametrisations for data and MC
  //
  
  //data type
  TString datatype="DATA";
  //in case of mc fRecoPass is per default 1
  if (fIsMC) {
    datatype="MC";
    fRecoPass=1;
  }
  
  //
  //set the PID splines
  //
  TString period=fLHCperiod;
  if (fArrPidResponseMaster){
    TObject *grAll=0x0;
    //for MC don't use period information
//     if (fIsMC) period="[A-Z0-9]*";
    //for MC use MC period information
    if (fIsMC) period=fMCperiodTPC;
//pattern for the default entry (valid for all particles)
    TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
    
    //loop over entries and filter them
    for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
      TObject *responseFunction=fArrPidResponseMaster->At(iresp);
      if (responseFunction==0x0) continue;
      TString responseName=responseFunction->GetName();
      
      if (!reg.MatchB(responseName)) continue;
      
      TObjArray *arr=reg.MatchS(responseName);
      TString particleName=arr->At(1)->GetName();
      delete arr;
      if (particleName.IsNull()) continue;
      if (particleName=="ALL") grAll=responseFunction;
      else {
        //find particle id
        for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
          TString particle=AliPID::ParticleName(ispec);
          particle.ToUpper();
          if ( particle == particleName ){
            //test if there is already a function set. If yes, cleanup
            TObject *old=const_cast<TObject*>(fPIDResponse->GetTPCResponse().GetResponseFunction((AliPID::EParticleType)ispec));
            if (old) delete old;
            fPIDResponse->GetTPCResponse().SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
            fPIDResponse->GetTPCResponse().SetUseDatabase(kTRUE);
            AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
            break;
          }
        }
      }
    }
    
    //set default response function to all particles which don't have a specific one
    if (grAll){
      for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
        if (!fPIDResponse->GetTPCResponse().GetResponseFunction((AliPID::EParticleType)ispec)){
          fPIDResponse->GetTPCResponse().SetResponseFunction((AliPID::EParticleType)ispec,grAll);
          AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
        }
      }
    }
  }

  //
  // Setup resolution parametrisation
  //

  //default
  fPIDResponse->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04);

  if (fRun>=122195){
    fPIDResponse->GetTPCResponse().SetSigma(2.30176e-02, 5.60422e+02);
  }
//   if ( fBeamType == "PBPB" ){
//     Double_t corrSigma=GetMultiplicityCorrectionSigma(GetTPCMultiplicityBin());
//     fPIDResponse->GetTPCResponse().SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
//   }
  
}

//______________________________________________________________________________
void AliAnalysisTaskPIDResponse::FillITSqa()
{
  AliVEvent *event=InputEvent();
  
  Int_t ntracks=event->GetNumberOfTracks();
  for(Int_t itrack = 0; itrack < ntracks; itrack++){
    AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
    ULong_t status=track->GetStatus();
    // not that nice. status bits not in virtual interface
    // ITS refit + ITS pid
    if (!( (status&0x0004==0x0004) && (status&0x0008==0x0008) )) return;
    Double_t mom=track->P();
    
    for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
      TH2 *h=(TH2*)fListQAits->At(ispecie);
      if (!h) continue;
      Double_t nSigma=fPIDResponse->NumberOfSigmasITS(track, (AliPID::EParticleType)ispecie);
      h->Fill(mom,nSigma);
    }
    
    TH2 *h=(TH2*)fListQAits->At(AliPID::kSPECIES);
    if (h) {
      Double_t sig=track->GetITSsignal();
      h->Fill(mom,sig);
    }
    
  }
}

//______________________________________________________________________________
void AliAnalysisTaskPIDResponse::FillTPCqa()
{
  AliVEvent *event=InputEvent();
  
  Int_t ntracks=event->GetNumberOfTracks();
  for(Int_t itrack = 0; itrack < ntracks; itrack++){
    AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
    
    //
    //basic track cuts
    //
    ULong_t status=track->GetStatus();
    // not that nice. status bits not in virtual interface
    // TPC refit + ITS refit + TPC pid
    if (!(status&0x0040==0x0040) || !(status&0x0004==0x0004) ||  !(status&0x0080==0x0080) ) return;
    
    Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
    Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
    if (track->GetTPCNclsF()>0) {
      ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
    }
    
    if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
    
    Double_t mom=track->GetTPCmomentum();
    
    for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
      TH2 *h=(TH2*)fListQAtpc->At(ispecie);
      if (!h) continue;
      Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
      h->Fill(mom,nSigma);
    }
    
    TH2 *h=(TH2*)fListQAtpc->At(AliPID::kSPECIES);
    if (h) {
      Double_t sig=track->GetTPCsignal();
      h->Fill(mom,sig);
    }
    
    TH2 *hPt=(TH2*)fListQAtpc->At(AliPID::kSPECIES+1);
    if (hPt) {
      Double_t sig=track->GetTPCsignal();
      hPt->Fill(sig);
    }
    
    TH2 *hMom=(TH2*)fListQAtpc->At(AliPID::kSPECIES+2);
    if (hMom) {
      Double_t sig=track->GetTPCmomentum();
      hMom->Fill(sig);
    }
    
  }
}

//______________________________________________________________________________
void AliAnalysisTaskPIDResponse::FillTOFqa()
{
  AliVEvent *event=InputEvent();
  
  Int_t ntracks=event->GetNumberOfTracks();
  for(Int_t itrack = 0; itrack < ntracks; itrack++){
    AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
    
    //
    //basic track cuts
    //
    ULong_t status=track->GetStatus();
    // not that nice. status bits not in virtual interface
    // TPC refit + ITS refit +
    // TOF out + TOFpid +
    // kTIME
    if (!(status&0x0040==0x0040) || !(status&0x0004==0x0004)
        || !(status&0x2000==0x2000) || !(status&0x8000==0x8000)
        || !(status&0x80000000==0x80000000) ) return;
    
    Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
    Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
    if (track->GetTPCNclsF()>0) {
      ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
    }
    
    if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
    
    
    Double_t mom=track->P();
    
    for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
      TH2 *h=(TH2*)fListQAtof->At(ispecie);
      if (!h) continue;
      Double_t nSigma=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
      h->Fill(mom,nSigma);
    }
    
    TH2 *h=(TH2*)fListQAtof->At(AliPID::kSPECIES);
    if (h) {
      Double_t sig=track->GetTOFsignal();
      h->Fill(mom,sig);
    }
    
  }
}

//______________________________________________________________________________
void AliAnalysisTaskPIDResponse::SetupTTSqa()
{
  //
  // Create the ITS qa objects
  //
  
  TVectorD *vX=MakeLogBinning(200,.1,30);
  
  for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
    TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_ITS_%s",AliPID::ParticleName(ispecie)),
                              Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
                              vX->GetNrows()-1,vX->GetMatrixArray(),
                              200,-10,10);
    fListQAits->Add(hNsigmaP);
  }
  
  
  TH2F *hSig = new TH2F("hSigP_ITS",
                        "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
                        vX->GetNrows()-1,vX->GetMatrixArray(),
                        300,0,300);
  
  fListQAits->Add(hSig);
  
}

//______________________________________________________________________________
void AliAnalysisTaskPIDResponse::SetupTPCqa()
{
  //
  // Create the TPC qa objects
  //
  
  TVectorD *vX=MakeLogBinning(200,.1,30);
  
  for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
    TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s",AliPID::ParticleName(ispecie)),
                              Form("TPC n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
                              vX->GetNrows()-1,vX->GetMatrixArray(),
                              200,-10,10);
    fListQAtpc->Add(hNsigmaP);
  }
  
  
  TH2F *hSig = new TH2F("hSigP_TPC",
                        "TPC signal vs. p;p [GeV]; TPC signal [arb. units]",
                        vX->GetNrows()-1,vX->GetMatrixArray(),
                        300,0,300);
  
  fListQAtpc->Add(hSig);
  
  TH1F *hPt = new TH1F("hPt","Pt_TPC",100,0,300);
  fListQAtpc->Add(hPt);
  
  TH1F *hMom = new TH1F("hMom","Mom_TPC",100,0,20);
  fListQAtpc->Add(hMom);
  
}

//______________________________________________________________________________
void AliAnalysisTaskPIDResponse::SetupTRDqa()
{
  //
  // Create the TRD qa objects
  //
  
}

//______________________________________________________________________________
void AliAnalysisTaskPIDResponse::SetupTOFqa()
{
  //
  // Create the TOF qa objects
  //
  
  TVectorD *vX=MakeLogBinning(200,.1,30);
  
  for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
    TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_%s",AliPID::ParticleName(ispecie)),
                              Form("TOF n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
                              vX->GetNrows()-1,vX->GetMatrixArray(),
                              200,-10,10);
    fListQAtof->Add(hNsigmaP);
  }
  
  
  TH2F *hSig = new TH2F("hSigP_TOF",
                        "TOF signal vs. p;p [GeV]; TOF signal [arb. units]",
                        vX->GetNrows()-1,vX->GetMatrixArray(),
                        300,0,300);
  
  fListQAtof->Add(hSig);
  
}

//______________________________________________________________________________
TVectorD* AliAnalysisTaskPIDResponse::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
{
  //
  // Make logarithmic binning
  // the user has to delete the array afterwards!!!
  //
  
  //check limits
  if (xmin<1e-20 || xmax<1e-20){
    AliError("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
    return MakeLinBinning(nbinsX, xmin, xmax);
  }
  if (xmax<xmin){
    Double_t tmp=xmin;
    xmin=xmax;
    xmax=tmp;
  }
  TVectorD *binLim=new TVectorD(nbinsX+1);
  Double_t first=xmin;
  Double_t last=xmax;
  Double_t expMax=TMath::Log(last/first);
  for (Int_t i=0; i<nbinsX+1; ++i){
    (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
  }
  return binLim;
}

//______________________________________________________________________________
TVectorD* AliAnalysisTaskPIDResponse::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
{
  //
  // Make linear binning
  // the user has to delete the array afterwards!!!
  //
  if (xmax<xmin){
    Double_t tmp=xmin;
    xmin=xmax;
    xmax=tmp;
  }
  TVectorD *binLim=new TVectorD(nbinsX+1);
  Double_t first=xmin;
  Double_t last=xmax;
  Double_t binWidth=(last-first)/nbinsX;
  for (Int_t i=0; i<nbinsX+1; ++i){
    (*binLim)[i]=first+binWidth*(Double_t)i;
  }
  return binLim;
}

//_____________________________________________________________________________
TVectorD* AliAnalysisTaskPIDResponse::MakeArbitraryBinning(const char* bins)
{
  //
  // Make arbitrary binning, bins separated by a ','
  //
  TString limits(bins);
  if (limits.IsNull()){
    AliError("Bin Limit string is empty, cannot add the variable");
    return 0x0;
  }
  
  TObjArray *arr=limits.Tokenize(",");
  Int_t nLimits=arr->GetEntries();
  if (nLimits<2){
    AliError("Need at leas 2 bin limits, cannot add the variable");
    delete arr;
    return 0x0;
  }
  
  TVectorD *binLimits=new TVectorD(nLimits);
  for (Int_t iLim=0; iLim<nLimits; ++iLim){
    (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
  }
  
  delete arr;
  return binLimits;
}
Commit	Line	Data
66ab8ab2	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, 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	/* $Id: AliAnalysisTaskPIDResponse.cxx 43811 2010-09-23 14:13:31Z wiechula $ */
	17	#include <TList.h>
	18	#include <TVectorD.h>
	19	#include <TObjArray.h>
	20	#include <TH2.h>
	21	#include <TFile.h>
	22	#include <TPRegexp.h>
	23	#include <TChain.h>
	24
	25	#include <AliAnalysisManager.h>
	26	#include <AliInputEventHandler.h>
	27	#include <AliVEventHandler.h>
	28	#include <AliVEvent.h>
	29	#include <AliVParticle.h>
	30	#include <AliVTrack.h>
	31	#include <AliLog.h>
	32	#include <AliPID.h>
	33	#include <AliPIDResponse.h>
	34	#include <AliITSPIDResponse.h>
	35	#include <AliTPCPIDResponse.h>
	36	#include <AliTRDPIDResponse.h>
	37	#include <AliTOFPIDResponse.h>
	38
	39	#include "AliAnalysisTaskPIDResponse.h"
	40
	41
	42	ClassImp(AliAnalysisTaskPIDResponse)
	43
	44	//______________________________________________________________________________
	45	AliAnalysisTaskPIDResponse::AliAnalysisTaskPIDResponse():
	46	AliAnalysisTaskSE(),
	47	fIsMC(kFALSE),
	48	fTOFTimeZeroTypeUser(-1),
	49	fTOFTimeZeroType(AliPIDResponse::kBest_T0),
	50	fTOFres(100.),
	51	fPIDResponse(0x0),
	52	fListQA(0x0),
	53	fListQAits(0x0),
	54	fListQAtpc(0x0),
	55	fListQAtrd(0x0),
	56	fListQAtof(0x0),
	57	fBeamType("PP"),
	58	fLHCperiod(),
	59	fMCperiodTPC(),
	60	fRecoPass(0),
	61	fRun(0),
	62	fOldRun(0),
	63	fArrPidResponseMaster(0x0)
	64	{
65	//
66	// Dummy constructor
67	//
68	}
69
70	//______________________________________________________________________________
71	AliAnalysisTaskPIDResponse::AliAnalysisTaskPIDResponse(const char* name):
72	AliAnalysisTaskSE(name),
73	fIsMC(kFALSE),
74	fTOFTimeZeroTypeUser(-1),
52c1c86c	75	fTOFTimeZeroType(AliPIDResponse::kBest_T0),
66ab8ab2	76	fTOFres(100.),
	77	fPIDResponse(0x0),
	78	fListQA(0x0),
	79	fListQAits(0x0),
	80	fListQAtpc(0x0),
	81	fListQAtrd(0x0),
	82	fListQAtof(0x0),
	83	fBeamType("PP"),
	84	fLHCperiod(),
	85	fMCperiodTPC(),
	86	fRecoPass(0),
	87	fRun(0),
	88	fOldRun(0),
	89	fArrPidResponseMaster(0x0)
	90	{
	91	//
	92	// Default constructor
	93	//
	94	DefineInput(0,TChain::Class());
	95	DefineOutput(1,TList::Class());
	96	}
	97
	98	//______________________________________________________________________________
	99	AliAnalysisTaskPIDResponse::~AliAnalysisTaskPIDResponse()
	100	{
	101	//
	102	// Destructor
	103	//
	104
	105	delete fArrPidResponseMaster;
	106	}
	107
	108	//______________________________________________________________________________
	109	void AliAnalysisTaskPIDResponse::ExecNewRun()
	110	{
	111	//
	112	// Things to Execute upon a new run
	113	//
	114
	115	SetRecoInfo();
	116
	117	SetITSParametrisation();
	118
	119	SetTPCPidResponseMaster();
	120	SetTPCParametrisation();
	121
	122	fPIDResponse->GetTOFResponse().SetTimeResolution(fTOFres);
	123	}
	124
	125	//______________________________________________________________________________
	126	void AliAnalysisTaskPIDResponse::UserCreateOutputObjects()
	127	{
	128	//
	129	// Create the output QA objects
	130	//
	131
	132	AliLog::SetClassDebugLevel("AliAnalysisTaskPIDResponse",10);
	133
	134	//input hander
	135	AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
	136	AliInputEventHandler inputHandler=dynamic_cast<AliInputEventHandler>(man->GetInputEventHandler());
	137	if (!inputHandler) AliFatal("Input handler needed");
	138
	139	//pid response object
140	inputHandler->CreatePIDResponse(fIsMC);
141	fPIDResponse=inputHandler->GetPIDResponse();
142
143	//
144	fListQA=new TList;
145	fListQAits=new TList;
146	fListQAits->SetName("ITS");
147	fListQAtpc=new TList;
148	fListQAtpc->SetName("TPC");
149	fListQAtrd=new TList;
150	fListQAtrd->SetName("TRD");
151	fListQAtof=new TList;
152	fListQAtof->SetName("TOF");
153
154	fListQA->Add(fListQAits);
155	fListQA->Add(fListQAtpc);
156	fListQA->Add(fListQAtrd);
157	fListQA->Add(fListQAtof);
158
159	SetupTTSqa();
160	SetupTPCqa();
161	SetupTRDqa();
162	SetupTOFqa();
163
164	PostData(1,fListQA);
165	}
166
167
168	//______________________________________________________________________________
169	void AliAnalysisTaskPIDResponse::UserExec(Option_t /option*/)
170	{
171	//
172	// Setup the PID response functions and fill the QA histograms
173	//
174
175	AliVEvent *event=InputEvent();
176	if (!event) return;
177	fRun=event->GetRunNumber();
178
179	if (fRun!=fOldRun){
180	ExecNewRun();
181	fOldRun=fRun;
182	}
183
184	Double_t timeZeroType=fTOFTimeZeroTypeUser;
185	if (timeZeroType<0) timeZeroType=fTOFTimeZeroType;
186	fPIDResponse->SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)timeZeroType);
187
188	FillITSqa();
189	FillTPCqa();
190	FillTOFqa();
191
192	PostData(1,fListQA);
193	}
194
195
196	//______________________________________________________________________________
197	void AliAnalysisTaskPIDResponse::SetRecoInfo()
198	{
199	//
200	// Set reconstruction information
201	//
202
203	//reset information
204	fRecoPass=0;
205	fLHCperiod="";
206	fMCperiodTPC="";
207
208	fBeamType="";
209
210	//Get the current file to check the reconstruction pass (UGLY, but not stored in ESD... )
211	AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
212	AliVEventHandler *inputHandler=mgr->GetInputEventHandler();
213	if (!inputHandler) return;
214
215	TTree tree= (TTree)inputHandler->GetTree();
216	TFile file= (TFile)tree->GetCurrentFile();
217
218	if (!file) {
219	AliError("Current file not found, cannot set reconstruction information");
220	return;
221	}
222
223	fBeamType="PP";
224
225	//find the period by run number (UGLY, but not stored in ESD and AOD... )
226	if (fRun>=114737&&fRun<=117223) { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1"; }
227	else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1"; }
228	else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
229	else if (fRun>=127719&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
230	else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
231	else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
232	else if (fRun>=136851&&fRun<=139517) { fLHCperiod="LHC10H"; fMCperiodTPC="LHC10H8"; fBeamType="PBPB"; }
233	else if (fRun>=139699) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
234
235	//find pass from file name (UGLY, but not stored in ESD... )
236	TString fileName(file->GetName());
237	if (fileName.Contains("/pass1")) {
238	fRecoPass=1;
239	} else if (fileName.Contains("/pass2")) {
240	fRecoPass=2;
241	} else if (fileName.Contains("/pass3")) {
242	fRecoPass=3;
243	}
244
245	}
246
247	//______________________________________________________________________________
248	void AliAnalysisTaskPIDResponse::SetITSParametrisation()
249	{
250	//
251	// Set the ITS parametrisation
252	//
253	}
254
255	//______________________________________________________________________________
256	void AliAnalysisTaskPIDResponse::SetTPCPidResponseMaster()
257	{
258	//
259	// Load the TPC pid response functions from the OADB
260	//
261
262	//reset the PID response functions
263	delete fArrPidResponseMaster;
264	fArrPidResponseMaster=0x0;
265
266	TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", AliAnalysisManager::GetOADBPath()));
267
268	TFile f(fileName.Data());
269	if (f.IsOpen() && !f.IsZombie()){
270	fArrPidResponseMaster=dynamic_cast<TObjArray*>(f.Get("TPCPIDResponse"));
271	fArrPidResponseMaster->SetOwner();
272	f.Close();
273	}
274
275	if (!fArrPidResponseMaster){
276	AliFatal("Could not retrieve the TPC pid response");
277	return;
278	}
279	}
280
281	//______________________________________________________________________________
282	void AliAnalysisTaskPIDResponse::SetTPCParametrisation()
283	{
284	//
285	// Change BB parametrisation for current run
286	//
287
288	if (fLHCperiod.IsNull()) {
289	AliFatal("No period set, not changing parametrisation");
290	return;
291	}
292
293	//
294	// Set default parametrisations for data and MC
295	//
296
297	//data type
298	TString datatype="DATA";
299	//in case of mc fRecoPass is per default 1
300	if (fIsMC) {
301	datatype="MC";
302	fRecoPass=1;
303	}
304
305	//
306	//set the PID splines
307	//
308	TString period=fLHCperiod;
309	if (fArrPidResponseMaster){
310	TObject *grAll=0x0;
311	//for MC don't use period information
312	// if (fIsMC) period="[A-Z0-9]*";
313	//for MC use MC period information
314	if (fIsMC) period=fMCperiodTPC;
315	//pattern for the default entry (valid for all particles)
316	TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
317
318	//loop over entries and filter them
319	for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
320	TObject *responseFunction=fArrPidResponseMaster->At(iresp);
321	if (responseFunction==0x0) continue;
322	TString responseName=responseFunction->GetName();
323
324	if (!reg.MatchB(responseName)) continue;
325
326	TObjArray *arr=reg.MatchS(responseName);
327	TString particleName=arr->At(1)->GetName();
328	delete arr;
329	if (particleName.IsNull()) continue;
330	if (particleName=="ALL") grAll=responseFunction;
331	else {
332	//find particle id
333	for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
334	TString particle=AliPID::ParticleName(ispec);
335	particle.ToUpper();
336	if ( particle == particleName ){
337	//test if there is already a function set. If yes, cleanup
338	TObject old=const_cast<TObject>(fPIDResponse->GetTPCResponse().GetResponseFunction((AliPID::EParticleType)ispec));
339	if (old) delete old;
340	fPIDResponse->GetTPCResponse().SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
341	fPIDResponse->GetTPCResponse().SetUseDatabase(kTRUE);
342	AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
343	break;
344	}
345	}
346	}
347	}
348
349	//set default response function to all particles which don't have a specific one
350	if (grAll){
351	for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
352	if (!fPIDResponse->GetTPCResponse().GetResponseFunction((AliPID::EParticleType)ispec)){
353	fPIDResponse->GetTPCResponse().SetResponseFunction((AliPID::EParticleType)ispec,grAll);
354	AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
355	}
356	}
357	}
358	}
359
360	//
361	// Setup resolution parametrisation
362	//
363
364	//default
365	fPIDResponse->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04);
366
367	if (fRun>=122195){
368	fPIDResponse->GetTPCResponse().SetSigma(2.30176e-02, 5.60422e+02);
369	}
370	// if ( fBeamType == "PBPB" ){
371	// Double_t corrSigma=GetMultiplicityCorrectionSigma(GetTPCMultiplicityBin());
372	// fPIDResponse->GetTPCResponse().SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
373	// }
374
375	}
376
377	//______________________________________________________________________________
378	void AliAnalysisTaskPIDResponse::FillITSqa()
379	{
380	AliVEvent *event=InputEvent();
381
382	Int_t ntracks=event->GetNumberOfTracks();
383	for(Int_t itrack = 0; itrack < ntracks; itrack++){
384	AliVTrack track=(AliVTrack)event->GetTrack(itrack);
385	ULong_t status=track->GetStatus();
386	// not that nice. status bits not in virtual interface
387	// ITS refit + ITS pid
388	if (!( (status&0x0004==0x0004) && (status&0x0008==0x0008) )) return;
389	Double_t mom=track->P();
390
391	for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
392	TH2 h=(TH2)fListQAits->At(ispecie);
393	if (!h) continue;
394	Double_t nSigma=fPIDResponse->NumberOfSigmasITS(track, (AliPID::EParticleType)ispecie);
395	h->Fill(mom,nSigma);
396	}
397
398	TH2 h=(TH2)fListQAits->At(AliPID::kSPECIES);
399	if (h) {
400	Double_t sig=track->GetITSsignal();
401	h->Fill(mom,sig);
402	}
403
404	}
405	}
406
407	//______________________________________________________________________________
408	void AliAnalysisTaskPIDResponse::FillTPCqa()
409	{
410	AliVEvent *event=InputEvent();
411
412	Int_t ntracks=event->GetNumberOfTracks();
413	for(Int_t itrack = 0; itrack < ntracks; itrack++){
414	AliVTrack track=(AliVTrack)event->GetTrack(itrack);
415
416	//
417	//basic track cuts
418	//
419	ULong_t status=track->GetStatus();
420	// not that nice. status bits not in virtual interface
421	// TPC refit + ITS refit + TPC pid
422	if (!(status&0x0040==0x0040) \|\| !(status&0x0004==0x0004) \|\| !(status&0x0080==0x0080) ) return;
423
424	Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
425	Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
426	if (track->GetTPCNclsF()>0) {
427	ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
428	}
429
430	if ( nCrossedRowsTPC<70 \|\| ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
431
432	Double_t mom=track->GetTPCmomentum();
433
434	for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
435	TH2 h=(TH2)fListQAtpc->At(ispecie);
436	if (!h) continue;
437	Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
438	h->Fill(mom,nSigma);
439	}
440
441	TH2 h=(TH2)fListQAtpc->At(AliPID::kSPECIES);
442	if (h) {
443	Double_t sig=track->GetTPCsignal();
444	h->Fill(mom,sig);
445	}
446
447	TH2 hPt=(TH2)fListQAtpc->At(AliPID::kSPECIES+1);
448	if (hPt) {
449	Double_t sig=track->GetTPCsignal();
450	hPt->Fill(sig);
451	}
452
453	TH2 hMom=(TH2)fListQAtpc->At(AliPID::kSPECIES+2);
454	if (hMom) {
455	Double_t sig=track->GetTPCmomentum();
456	hMom->Fill(sig);
457	}
458
459	}
460	}
461
462	//______________________________________________________________________________
463	void AliAnalysisTaskPIDResponse::FillTOFqa()
464	{
465	AliVEvent *event=InputEvent();
466
467	Int_t ntracks=event->GetNumberOfTracks();
468	for(Int_t itrack = 0; itrack < ntracks; itrack++){
469	AliVTrack track=(AliVTrack)event->GetTrack(itrack);
470
471	//
472	//basic track cuts
473	//
474	ULong_t status=track->GetStatus();
475	// not that nice. status bits not in virtual interface
476	// TPC refit + ITS refit +
477	// TOF out + TOFpid +
478	// kTIME
479	if (!(status&0x0040==0x0040) \|\| !(status&0x0004==0x0004)
480	\|\| !(status&0x2000==0x2000) \|\| !(status&0x8000==0x8000)
481	\|\| !(status&0x80000000==0x80000000) ) return;
482
483	Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
484	Float_t ratioCrossedRowsOverFindableClustersTPC = 1.0;
485	if (track->GetTPCNclsF()>0) {
486	ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
487	}
488
489	if ( nCrossedRowsTPC<70 \|\| ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
490
491
492	Double_t mom=track->P();
493
494	for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
495	TH2 h=(TH2)fListQAtof->At(ispecie);
496	if (!h) continue;
497	Double_t nSigma=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
498	h->Fill(mom,nSigma);
499	}
500
501	TH2 h=(TH2)fListQAtof->At(AliPID::kSPECIES);
502	if (h) {
503	Double_t sig=track->GetTOFsignal();
504	h->Fill(mom,sig);
505	}
506
507	}
508	}
509
510	//______________________________________________________________________________
511	void AliAnalysisTaskPIDResponse::SetupTTSqa()
512	{
513	//
514	// Create the ITS qa objects
515	//
516
517	TVectorD *vX=MakeLogBinning(200,.1,30);
518
519	for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
520	TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_ITS_%s",AliPID::ParticleName(ispecie)),
521	Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
522	vX->GetNrows()-1,vX->GetMatrixArray(),
523	200,-10,10);
524	fListQAits->Add(hNsigmaP);
525	}
526
527
528	TH2F *hSig = new TH2F("hSigP_ITS",
529	"ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
530	vX->GetNrows()-1,vX->GetMatrixArray(),
531	300,0,300);
532
533	fListQAits->Add(hSig);
534
535	}
536
537	//______________________________________________________________________________
538	void AliAnalysisTaskPIDResponse::SetupTPCqa()
539	{
540	//
541	// Create the TPC qa objects
542	//
543
544	TVectorD *vX=MakeLogBinning(200,.1,30);
545
546	for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
547	TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s",AliPID::ParticleName(ispecie)),
548	Form("TPC n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
549	vX->GetNrows()-1,vX->GetMatrixArray(),
550	200,-10,10);
551	fListQAtpc->Add(hNsigmaP);
552	}
553
554
555	TH2F *hSig = new TH2F("hSigP_TPC",
556	"TPC signal vs. p;p [GeV]; TPC signal [arb. units]",
557	vX->GetNrows()-1,vX->GetMatrixArray(),
558	300,0,300);
559
560	fListQAtpc->Add(hSig);
561
562	TH1F *hPt = new TH1F("hPt","Pt_TPC",100,0,300);
563	fListQAtpc->Add(hPt);
564
565	TH1F *hMom = new TH1F("hMom","Mom_TPC",100,0,20);
566	fListQAtpc->Add(hMom);
567
568	}
569
570	//______________________________________________________________________________
571	void AliAnalysisTaskPIDResponse::SetupTRDqa()
572	{
573	//
574	// Create the TRD qa objects
575	//
576
577	}
578
579	//______________________________________________________________________________
580	void AliAnalysisTaskPIDResponse::SetupTOFqa()
581	{
582	//
583	// Create the TOF qa objects
584	//
585
586	TVectorD *vX=MakeLogBinning(200,.1,30);
587
588	for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
589	TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_%s",AliPID::ParticleName(ispecie)),
590	Form("TOF n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
591	vX->GetNrows()-1,vX->GetMatrixArray(),
592	200,-10,10);
593	fListQAtof->Add(hNsigmaP);
594	}
595
596
597	TH2F *hSig = new TH2F("hSigP_TOF",
598	"TOF signal vs. p;p [GeV]; TOF signal [arb. units]",
599	vX->GetNrows()-1,vX->GetMatrixArray(),
600	300,0,300);
601
602	fListQAtof->Add(hSig);
603
604	}
605
606	//______________________________________________________________________________
607	TVectorD* AliAnalysisTaskPIDResponse::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
608	{
609	//
610	// Make logarithmic binning
611	// the user has to delete the array afterwards!!!
612	//
613
614	//check limits
615	if (xmin<1e-20 \|\| xmax<1e-20){
616	AliError("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
617	return MakeLinBinning(nbinsX, xmin, xmax);
618	}
619	if (xmax<xmin){
620	Double_t tmp=xmin;
621	xmin=xmax;
622	xmax=tmp;
623	}
624	TVectorD *binLim=new TVectorD(nbinsX+1);
625	Double_t first=xmin;
626	Double_t last=xmax;
627	Double_t expMax=TMath::Log(last/first);
628	for (Int_t i=0; i<nbinsX+1; ++i){
629	(binLim)[i]=firstTMath::Exp(expMax/nbinsX*(Double_t)i);
630	}
631	return binLim;
632	}
633
634	//______________________________________________________________________________
635	TVectorD* AliAnalysisTaskPIDResponse::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
636	{
637	//
638	// Make linear binning
639	// the user has to delete the array afterwards!!!
640	//
641	if (xmax<xmin){
642	Double_t tmp=xmin;
643	xmin=xmax;
644	xmax=tmp;
645	}
646	TVectorD *binLim=new TVectorD(nbinsX+1);
647	Double_t first=xmin;
648	Double_t last=xmax;
649	Double_t binWidth=(last-first)/nbinsX;
650	for (Int_t i=0; i<nbinsX+1; ++i){
651	(binLim)[i]=first+binWidth(Double_t)i;
652	}
653	return binLim;
654	}
655
656	//_____________________________________________________________________________
657	TVectorD* AliAnalysisTaskPIDResponse::MakeArbitraryBinning(const char* bins)
658	{
659	//
660	// Make arbitrary binning, bins separated by a ','
661	//
662	TString limits(bins);
663	if (limits.IsNull()){
664	AliError("Bin Limit string is empty, cannot add the variable");
665	return 0x0;
666	}
667
668	TObjArray *arr=limits.Tokenize(",");
669	Int_t nLimits=arr->GetEntries();
670	if (nLimits<2){
671	AliError("Need at leas 2 bin limits, cannot add the variable");
672	delete arr;
673	return 0x0;
674	}
675
676	TVectorD *binLimits=new TVectorD(nLimits);
677	for (Int_t iLim=0; iLim<nLimits; ++iLim){
678	(binLimits)[iLim]=(static_cast<TObjString>(arr->At(iLim)))->GetString().Atof();
679	}
680
681	delete arr;
682	return binLimits;
683	}