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   if (!fPIDResponse) AliFatal("PIDResponse object was not created");
 143
 144   //
 145   fListQA=new TList;
 146   fListQA->SetOwner();
 147
 148   fListQAits=new TList;
 149   fListQAits->SetOwner();
 150   fListQAits->SetName("ITS");
 151
 152   fListQAtpc=new TList;
 153   fListQAtpc->SetOwner();
 154   fListQAtpc->SetName("TPC");
 155
 156   fListQAtrd=new TList;
 157   fListQAtrd->SetOwner();
 158   fListQAtrd->SetName("TRD");
 159
 160   fListQAtof=new TList;
 161   fListQAtof->SetOwner();
 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"));
 281     f.Close();
 282   }
 283
 284   if (!fArrPidResponseMaster){
 285     AliFatal("Could not retrieve the TPC pid response");
 286     return;
 287   }
 288   fArrPidResponseMaster->SetOwner();
 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]=first*TMath::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 }