ANALYSIS/AliAnalysisTaskPIDResponse.cxx

   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),
  75 fTOFTimeZeroType(AliPIDResponse::kBest_T0),
  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]=first*TMath::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 }