[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();
  if (!fPIDResponse) AliFatal("PIDResponse object was not created");
  
  //
  fListQA=new TList;
  fListQA->SetOwner();
  
  fListQAits=new TList;
  fListQAits->SetOwner();
  fListQAits->SetName("ITS");
  
  fListQAtpc=new TList;
  fListQAtpc->SetOwner();
  fListQAtpc->SetName("TPC");
  
  fListQAtrd=new TList;
  fListQAtrd->SetOwner();
  fListQAtrd->SetName("TRD");
  
  fListQAtof=new TList;
  fListQAtof->SetOwner();
  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"));
    f.Close();
  }

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

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