ANALYSIS/AliAnalysisTaskPIDqa.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: AliAnalysisTaskPIDqa.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 #include <TF1.h>
  25 #include <TSpline.h>
  26
  27 #include <AliAnalysisManager.h>
  28 #include <AliInputEventHandler.h>
  29 #include <AliVEventHandler.h>
  30 #include <AliVEvent.h>
  31 #include <AliVParticle.h>
  32 #include <AliVTrack.h>
  33 #include <AliLog.h>
  34 #include <AliPID.h>
  35 #include <AliPIDResponse.h>
  36 #include <AliITSPIDResponse.h>
  37 #include <AliTPCPIDResponse.h>
  38 #include <AliTRDPIDResponse.h>
  39 #include <AliTOFPIDResponse.h>
  40
  41 #include <AliESDEvent.h>
  42
  43 #include "AliAnalysisTaskPIDqa.h"
  44
  45
  46 ClassImp(AliAnalysisTaskPIDqa)
  47
  48 //______________________________________________________________________________
  49 AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa():
  50 AliAnalysisTaskSE(),
  51 fPIDResponse(0x0),
  52 fListQA(0x0),
  53 fListQAits(0x0),
  54 fListQAitsSA(0x0),
  55 fListQAitsPureSA(0x0),
  56 fListQAtpc(0x0),
  57 fListQAtrd(0x0),
  58 fListQAtof(0x0),
  59 fListQAemcal(0x0),
  60 fListQAtpctof(0x0)
  61 {
  62   //
  63   // Dummy constructor
  64   //
  65 }
  66
  67 //______________________________________________________________________________
  68 AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa(const char* name):
  69 AliAnalysisTaskSE(name),
  70 fPIDResponse(0x0),
  71 fListQA(0x0),
  72 fListQAits(0x0),
  73 fListQAitsSA(0x0),
  74 fListQAitsPureSA(0x0),
  75 fListQAtpc(0x0),
  76 fListQAtrd(0x0),
  77 fListQAtof(0x0),
  78 fListQAemcal(0x0),
  79 fListQAtpctof(0x0)
  80 {
  81   //
  82   // Default constructor
  83   //
  84   DefineInput(0,TChain::Class());
  85   DefineOutput(1,TList::Class());
  86 }
  87
  88 //______________________________________________________________________________
  89 AliAnalysisTaskPIDqa::~AliAnalysisTaskPIDqa()
  90 {
  91   //
  92   // Destructor
  93   //
  94
  95 }
  96
  97 //______________________________________________________________________________
  98 void AliAnalysisTaskPIDqa::UserCreateOutputObjects()
  99 {
 100   //
 101   // Create the output QA objects
 102   //
 103
 104   AliLog::SetClassDebugLevel("AliAnalysisTaskPIDqa",10);
 105
 106   //input hander
 107   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 108   AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
 109   if (!inputHandler) AliFatal("Input handler needed");
 110
 111   //pid response object
 112   fPIDResponse=inputHandler->GetPIDResponse();
 113   if (!fPIDResponse) AliError("PIDResponse object was not created");
 114
 115   //
 116   fListQA=new TList;
 117   fListQA->SetOwner();
 118
 119   fListQAits=new TList;
 120   fListQAits->SetOwner();
 121   fListQAits->SetName("ITS");
 122
 123   fListQAitsSA=new TList;
 124   fListQAitsSA->SetOwner();
 125   fListQAitsSA->SetName("ITS");
 126
 127   fListQAitsPureSA=new TList;
 128   fListQAitsPureSA->SetOwner();
 129   fListQAitsPureSA->SetName("ITS");
 130
 131   fListQAtpc=new TList;
 132   fListQAtpc->SetOwner();
 133   fListQAtpc->SetName("TPC");
 134
 135   fListQAtrd=new TList;
 136   fListQAtrd->SetOwner();
 137   fListQAtrd->SetName("TRD");
 138
 139   fListQAtof=new TList;
 140   fListQAtof->SetOwner();
 141   fListQAtof->SetName("TOF");
 142
 143   fListQAemcal=new TList;
 144   fListQAemcal->SetOwner();
 145   fListQAemcal->SetName("EMCAL");
 146
 147   fListQAtpctof=new TList;
 148   fListQAtpctof->SetOwner();
 149   fListQAtpctof->SetName("TPC_TOF");
 150
 151   fListQA->Add(fListQAits);
 152   fListQA->Add(fListQAitsSA);
 153   fListQA->Add(fListQAitsPureSA);
 154   fListQA->Add(fListQAtpc);
 155   fListQA->Add(fListQAtrd);
 156   fListQA->Add(fListQAtof);
 157   fListQA->Add(fListQAemcal);
 158   fListQA->Add(fListQAtpctof);
 159
 160   SetupITSqa();
 161   SetupTPCqa();
 162   SetupTRDqa();
 163   SetupTOFqa();
 164   SetupEMCALqa();
 165   SetupTPCTOFqa();
 166
 167   PostData(1,fListQA);
 168 }
 169
 170
 171 //______________________________________________________________________________
 172 void AliAnalysisTaskPIDqa::UserExec(Option_t */*option*/)
 173 {
 174   //
 175   // Setup the PID response functions and fill the QA histograms
 176   //
 177
 178   AliVEvent *event=InputEvent();
 179   if (!event||!fPIDResponse) return;
 180
 181
 182   FillITSqa();
 183   FillTPCqa();
 184   FillTRDqa();
 185   FillTOFqa();
 186   FillEMCALqa();
 187   FillTPCTOFqa();
 188
 189   PostData(1,fListQA);
 190 }
 191
 192 //______________________________________________________________________________
 193 void AliAnalysisTaskPIDqa::FillITSqa()
 194 {
 195   //
 196   // Fill PID qa histograms for the ITS
 197   //
 198
 199   AliVEvent *event=InputEvent();
 200
 201   Int_t ntracks=event->GetNumberOfTracks();
 202   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 203     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 204     ULong_t status=track->GetStatus();
 205     // not that nice. status bits not in virtual interface
 206     // ITS refit + ITS pid selection
 207     if (!( ( (status & AliVTrack::kITSrefit)==AliVTrack::kITSrefit ) ||
 208            ! ( (status & AliVTrack::kITSpid  )==AliVTrack::kITSpid   ) )) continue;
 209     Double_t mom=track->P();
 210
 211     TList *theList = 0x0;
 212     if(( (status & AliVTrack::kTPCin)==AliVTrack::kTPCin )){
 213       //ITS+TPC tracks
 214       theList=fListQAits;
 215     }else{
 216       if(!( (status & AliVTrack::kITSpureSA)==AliVTrack::kITSpureSA )){
 217         //ITS Standalone tracks
 218         theList=fListQAitsSA;
 219       }else{
 220         //ITS Pure Standalone tracks
 221         theList=fListQAitsPureSA;
 222       }
 223     }
 224
 225
 226     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 227       TH2 *h=(TH2*)theList->At(ispecie);
 228       if (!h) continue;
 229       Double_t nSigma=fPIDResponse->NumberOfSigmasITS(track, (AliPID::EParticleType)ispecie);
 230       h->Fill(mom,nSigma);
 231     }
 232     TH2 *h=(TH2*)theList->At(AliPID::kSPECIES);
 233     if (h) {
 234       Double_t sig=track->GetITSsignal();
 235       h->Fill(mom,sig);
 236     }
 237   }
 238 }
 239
 240 //______________________________________________________________________________
 241 void AliAnalysisTaskPIDqa::FillTPCqa()
 242 {
 243   //
 244   // Fill PID qa histograms for the TPC
 245   //
 246
 247   AliVEvent *event=InputEvent();
 248
 249   Int_t ntracks=event->GetNumberOfTracks();
 250   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 251     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 252
 253     //
 254     //basic track cuts
 255     //
 256     ULong_t status=track->GetStatus();
 257     // not that nice. status bits not in virtual interface
 258     // TPC refit + ITS refit + TPC pid
 259     if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 260         !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
 261         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  ) ) continue;
 262
 263     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 264     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 265     if (track->GetTPCNclsF()>0) {
 266       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 267     }
 268
 269     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 270
 271     Double_t mom=track->GetTPCmomentum();
 272
 273     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 274       TH2 *h=(TH2*)fListQAtpc->At(ispecie);
 275       if (!h) continue;
 276       Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
 277       h->Fill(mom,nSigma);
 278     }
 279
 280     TH2 *h=(TH2*)fListQAtpc->At(AliPID::kSPECIES);
 281     if (h) {
 282       Double_t sig=track->GetTPCsignal();
 283       h->Fill(mom,sig);
 284     }
 285   }
 286 }
 287
 288 //______________________________________________________________________________
 289 void AliAnalysisTaskPIDqa::FillTRDqa()
 290 {
 291   //
 292   // Fill PID qa histograms for the TRD
 293   //
 294
 295 }
 296
 297 //______________________________________________________________________________
 298 void AliAnalysisTaskPIDqa::FillTOFqa()
 299 {
 300   //
 301   // Fill PID qa histograms for the TOF
 302   //
 303
 304   AliVEvent *event=InputEvent();
 305
 306
 307   Int_t ntracks=event->GetNumberOfTracks();
 308   Int_t tracksAtTof = 0;
 309   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 310     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 311
 312     //
 313     //basic track cuts
 314     //
 315     ULong_t status=track->GetStatus();
 316     // TPC refit + ITS refit +
 317     // TOF out + TOFpid +
 318     // kTIME
 319     // (we don't use kTOFmismatch because it depends on TPC....)
 320     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 321         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
 322         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
 323         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
 324         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
 325
 326     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 327     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 328     if (track->GetTPCNclsF()>0) {
 329       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 330     }
 331
 332     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 333
 334     tracksAtTof++;
 335
 336     Double_t mom=track->P();
 337
 338     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 339       TH2 *h=(TH2*)fListQAtof->At(ispecie);
 340       if (!h) continue;
 341       Double_t nSigma=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
 342       h->Fill(mom,nSigma);
 343     }
 344
 345     TH2 *h=(TH2*)fListQAtof->FindObject("hSigP_TOF");
 346     if (h) {
 347       Double_t sig=track->GetTOFsignal();
 348       h->Fill(mom,sig);
 349     }
 350
 351     Double_t mask = (Double_t)fPIDResponse->GetTOFResponse().GetStartTimeMask(mom) + 0.5;
 352     ((TH1F*)fListQAtof->FindObject("hStartTimeMask_TOF"))->Fill(mask);
 353
 354     Double_t res = (Double_t)fPIDResponse->GetTOFResponse().GetStartTimeRes(mom);
 355     ((TH1F*)fListQAtof->FindObject("hStartTimeRes_TOF"))->Fill(res);
 356
 357     AliESDEvent *esd = dynamic_cast<AliESDEvent *>(event);
 358     if (esd) {
 359       Double_t startTime = esd->GetT0TOF(0);
 360       if (startTime < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeAC_T0"))->Fill(startTime);
 361       else {
 362         startTime = esd->GetT0TOF(1);
 363         if (startTime < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeA_T0"))->Fill(startTime);
 364         startTime = esd->GetT0TOF(2);
 365         if (startTime < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeC_T0"))->Fill(startTime);
 366       }
 367     }
 368   }
 369   if (tracksAtTof > 0) {
 370     ((TH1F* )fListQAtof->FindObject("hnTracksAt_TOF"))->Fill(tracksAtTof);
 371     Int_t mask = fPIDResponse->GetTOFResponse().GetStartTimeMask(5.);
 372     if (mask & 0x1) ((TH1F*)fListQAtof->FindObject("hT0MakerEff"))->Fill(tracksAtTof);
 373   }
 374
 375 }
 376
 377 //______________________________________________________________________________
 378 void AliAnalysisTaskPIDqa::FillEMCALqa()
 379 {
 380   //
 381   // Fill PID qa histograms for the EMCAL
 382   //
 383
 384   AliVEvent *event=InputEvent();
 385
 386   Int_t ntracks=event->GetNumberOfTracks();
 387   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 388     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 389
 390     //
 391     //basic track cuts
 392     //
 393     ULong_t status=track->GetStatus();
 394     // not that nice. status bits not in virtual interface
 395     // TPC refit + ITS refit +
 396     // TOF out + TOFpid +
 397     // kTIME
 398     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
 399
 400     Double_t pt=track->Pt();
 401
 402     //EMCAL nSigma (only for electrons at the moment)
 403     TH2 *h=(TH2*)fListQAemcal->At(0);
 404     if (!h) continue;
 405     Double_t nSigma=fPIDResponse->NumberOfSigmasEMCAL(track, (AliPID::EParticleType)0);
 406     h->Fill(pt,nSigma);
 407
 408     //EMCAL signal (E/p vs. pT)
 409     h=(TH2*)fListQAemcal->At(1);
 410     if (h) {
 411
 412       Int_t nMatchClus = track->GetEMCALcluster();
 413       Double_t mom     = track->P();
 414       Double_t eop     = -1.;
 415
 416       if(nMatchClus > -1){
 417
 418         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
 419
 420         if(matchedClus){
 421
 422           // matched cluster is EMCAL
 423           if(matchedClus->IsEMCAL()){
 424
 425             Double_t fClsE       = matchedClus->E();
 426             eop                  = fClsE/mom;
 427
 428             h->Fill(pt,eop);
 429
 430           }
 431         }
 432       }
 433       else{
 434         Printf("status status = AliVTrack::kEMCALmatch, BUT no matched cluster!");
 435       }
 436     }
 437   }
 438 }
 439
 440 //______________________________________________________________________________
 441 void AliAnalysisTaskPIDqa::FillTPCTOFqa()
 442 {
 443   //
 444   // Fill PID qa histograms for the TOF
 445   //   Here also the TPC histograms after TOF selection are filled
 446   //
 447
 448   AliVEvent *event=InputEvent();
 449
 450   Int_t ntracks=event->GetNumberOfTracks();
 451   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 452     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 453
 454     //
 455     //basic track cuts
 456     //
 457     ULong_t status=track->GetStatus();
 458     // not that nice. status bits not in virtual interface
 459     // TPC refit + ITS refit +
 460     // TOF out + TOFpid +
 461     // kTIME
 462     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 463         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
 464         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid ) ||
 465         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
 466         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
 467         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
 468
 469     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 470     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 471     if (track->GetTPCNclsF()>0) {
 472       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 473     }
 474
 475     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 476
 477
 478     Double_t mom=track->P();
 479     Double_t momTPC=track->GetTPCmomentum();
 480
 481     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 482       //TOF nSigma
 483       Double_t nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
 484       Double_t nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
 485
 486       //TPC after TOF cut
 487       TH2 *h=(TH2*)fListQAtpctof->At(ispecie);
 488       if (h && TMath::Abs(nSigmaTOF)<3.) h->Fill(momTPC,nSigmaTPC);
 489
 490       //TOF after TPC cut
 491       h=(TH2*)fListQAtpctof->At(ispecie+AliPID::kSPECIES);
 492       if (h && TMath::Abs(nSigmaTPC)<3.) h->Fill(mom,nSigmaTOF);
 493
 494       //EMCAL after TOF and TPC cut
 495       h=(TH2*)fListQAtpctof->At(ispecie+2*AliPID::kSPECIES);
 496       if (h && TMath::Abs(nSigmaTOF)<3. && TMath::Abs(nSigmaTPC)<3. ){
 497
 498         Int_t nMatchClus = track->GetEMCALcluster();
 499         Double_t pt      = track->Pt();
 500         Double_t eop     = -1.;
 501
 502         if(nMatchClus > -1){
 503
 504           AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
 505
 506           if(matchedClus){
 507
 508             // matched cluster is EMCAL
 509             if(matchedClus->IsEMCAL()){
 510
 511               Double_t fClsE       = matchedClus->E();
 512               eop                  = fClsE/mom;
 513
 514               h->Fill(pt,eop);
 515
 516
 517             }
 518           }
 519         }
 520       }
 521     }
 522   }
 523 }
 524
 525 //______________________________________________________________________________
 526 void AliAnalysisTaskPIDqa::SetupITSqa()
 527 {
 528   //
 529   // Create the ITS qa objects
 530   //
 531
 532   TVectorD *vX=MakeLogBinning(200,.1,30);
 533
 534   //ITS+TPC tracks
 535   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 536     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_ITS_%s",AliPID::ParticleName(ispecie)),
 537                               Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 538                               vX->GetNrows()-1,vX->GetMatrixArray(),
 539                               200,-10,10);
 540     fListQAits->Add(hNsigmaP);
 541   }
 542   TH2F *hSig = new TH2F("hSigP_ITS",
 543                         "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
 544                         vX->GetNrows()-1,vX->GetMatrixArray(),
 545                         300,0,300);
 546   fListQAits->Add(hSig);
 547
 548   //ITS Standalone tracks
 549   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 550     TH2F *hNsigmaPSA = new TH2F(Form("hNsigmaP_ITSSA_%s",AliPID::ParticleName(ispecie)),
 551                                 Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 552                                 vX->GetNrows()-1,vX->GetMatrixArray(),
 553                                 200,-10,10);
 554     fListQAitsSA->Add(hNsigmaPSA);
 555   }
 556   TH2F *hSigSA = new TH2F("hSigP_ITSSA",
 557                           "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
 558                           vX->GetNrows()-1,vX->GetMatrixArray(),
 559                           300,0,300);
 560   fListQAitsSA->Add(hSigSA);
 561
 562   //ITS Pure Standalone tracks
 563   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 564     TH2F *hNsigmaPPureSA = new TH2F(Form("hNsigmaP_ITSPureSA_%s",AliPID::ParticleName(ispecie)),
 565                                     Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 566                                     vX->GetNrows()-1,vX->GetMatrixArray(),
 567                                     200,-10,10);
 568     fListQAitsPureSA->Add(hNsigmaPPureSA);
 569   }
 570   TH2F *hSigPureSA = new TH2F("hSigP_ITSPureSA",
 571                               "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
 572                               vX->GetNrows()-1,vX->GetMatrixArray(),
 573                               300,0,300);
 574   fListQAitsPureSA->Add(hSigPureSA);
 575
 576   delete vX;
 577 }
 578
 579 //______________________________________________________________________________
 580 void AliAnalysisTaskPIDqa::SetupTPCqa()
 581 {
 582   //
 583   // Create the TPC 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_TPC_%s",AliPID::ParticleName(ispecie)),
 590                               Form("TPC n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 591                               vX->GetNrows()-1,vX->GetMatrixArray(),
 592                               200,-10,10);
 593     fListQAtpc->Add(hNsigmaP);
 594   }
 595
 596
 597   TH2F *hSig = new TH2F("hSigP_TPC",
 598                         "TPC signal vs. p;p [GeV]; TPC signal [arb. units]",
 599                         vX->GetNrows()-1,vX->GetMatrixArray(),
 600                         300,0,300);
 601   fListQAtpc->Add(hSig);
 602
 603   delete vX;
 604 }
 605
 606 //______________________________________________________________________________
 607 void AliAnalysisTaskPIDqa::SetupTRDqa()
 608 {
 609   //
 610   // Create the TRD qa objects
 611   //
 612
 613 }
 614
 615 //______________________________________________________________________________
 616 void AliAnalysisTaskPIDqa::SetupTOFqa()
 617 {
 618   //
 619   // Create the TOF qa objects
 620   //
 621
 622   TVectorD *vX=MakeLogBinning(200,.1,30);
 623
 624   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 625     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_%s",AliPID::ParticleName(ispecie)),
 626                               Form("TOF n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 627                               vX->GetNrows()-1,vX->GetMatrixArray(),
 628                               200,-10,10);
 629     fListQAtof->Add(hNsigmaP);
 630     // to be added: t-texp without StartTime subtraction
 631   }
 632
 633
 634   TH2F *hSig = new TH2F("hSigP_TOF",
 635                         "TOF signal vs. p;p [GeV]; TOF signal [arb. units]",
 636                         vX->GetNrows()-1,vX->GetMatrixArray(),
 637                         300,0,300);
 638
 639   delete vX;
 640
 641   fListQAtof->Add(hSig);
 642
 643   TH1F *hStartTimeMask_TOF = new TH1F("hStartTimeMask_TOF","StartTime mask",8,0,8);
 644   fListQAtof->Add(hStartTimeMask_TOF);
 645   TH1F *hStartTimeRes_TOF = new TH1F("hStartTimeRes_TOF","StartTime resolution [ps]",100,0,500);
 646   fListQAtof->Add(hStartTimeRes_TOF);
 647
 648   TH1F *hnTracksAt_TOF = new TH1F("hnTracksAt_TOF","Matched tracks at TOF",20,0,20);
 649   fListQAtof->Add(hnTracksAt_TOF);
 650   TH1F *hT0MakerEff = new TH1F("hT0MakerEff","Events with T0-TOF vs nTracks",20,0,20);
 651   fListQAtof->Add(hT0MakerEff);
 652
 653   // this in principle should stay on a T0 PID QA, but are just the data prepared for TOF use
 654   TH1F *hStartTimeA_T0 = new TH1F("hStartTimeA_T0","StartTime from T0A [ps]",1000,-1000,1000);
 655   fListQAtof->Add(hStartTimeA_T0);
 656   TH1F *hStartTimeC_T0 = new TH1F("hStartTimeC_T0","StartTime from T0C [ps]",1000,-1000,1000);
 657   fListQAtof->Add(hStartTimeC_T0);
 658   TH1F *hStartTimeAC_T0 = new TH1F("hStartTimeAC_T0","StartTime from T0AC [ps]",1000,-1000,1000);;
 659   fListQAtof->Add(hStartTimeAC_T0);
 660 }
 661
 662 //______________________________________________________________________________
 663 void AliAnalysisTaskPIDqa::SetupEMCALqa()
 664 {
 665   //
 666   // Create the EMCAL qa objects
 667   //
 668
 669   TVectorD *vX=MakeLogBinning(200,.1,30);
 670
 671   TH2F *hNsigmaPt = new TH2F(Form("hNsigmaPt_EMCAL_%s",AliPID::ParticleName(0)),
 672                              Form("EMCAL n#sigma %s vs. p_{T};p_{T} [GeV]; n#sigma",AliPID::ParticleName(0)),
 673                              vX->GetNrows()-1,vX->GetMatrixArray(),
 674                              200,-10,10);
 675   fListQAemcal->Add(hNsigmaPt);
 676
 677   TH2F *hSigPt = new TH2F("hSigPt_EMCAL",
 678                         "EMCAL signal (E/p) vs. p_{T};p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
 679                         vX->GetNrows()-1,vX->GetMatrixArray(),
 680                         200,0,2);
 681   fListQAemcal->Add(hSigPt);
 682
 683   delete vX;
 684 }
 685
 686 //______________________________________________________________________________
 687 void AliAnalysisTaskPIDqa::SetupTPCTOFqa()
 688 {
 689   //
 690   // Create the qa objects for TPC + TOF combination
 691   //
 692
 693   TVectorD *vX=MakeLogBinning(200,.1,30);
 694
 695   //TPC signals after TOF cut
 696   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 697     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_TOF_%s",AliPID::ParticleName(ispecie)),
 698                               Form("TPC n#sigma %s vs. p (after TOF 3#sigma cut);p_{TPC} [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 699                               vX->GetNrows()-1,vX->GetMatrixArray(),
 700                               200,-10,10);
 701     fListQAtpctof->Add(hNsigmaP);
 702   }
 703
 704   //TOF signals after TPC cut
 705   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 706     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
 707                               Form("TOF n#sigma %s vs. p (after TPC n#sigma cut);p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 708                               vX->GetNrows()-1,vX->GetMatrixArray(),
 709                               200,-10,10);
 710     fListQAtpctof->Add(hNsigmaP);
 711   }
 712
 713   //EMCAL signal after TOF and TPC cut
 714   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 715     TH2F *heopPt = new TH2F(Form("heopPt_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
 716                             Form("EMCAL signal (E/p) %s vs. p_{T};p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",AliPID::ParticleName(ispecie)),
 717                             vX->GetNrows()-1,vX->GetMatrixArray(),
 718                             200,0,2);
 719     fListQAtpctof->Add(heopPt);
 720   }
 721
 722   delete vX;
 723 }
 724
 725 //______________________________________________________________________________
 726 TVectorD* AliAnalysisTaskPIDqa::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
 727 {
 728   //
 729   // Make logarithmic binning
 730   // the user has to delete the array afterwards!!!
 731   //
 732
 733   //check limits
 734   if (xmin<1e-20 || xmax<1e-20){
 735     AliError("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
 736     return MakeLinBinning(nbinsX, xmin, xmax);
 737   }
 738   if (xmax<xmin){
 739     Double_t tmp=xmin;
 740     xmin=xmax;
 741     xmax=tmp;
 742   }
 743   TVectorD *binLim=new TVectorD(nbinsX+1);
 744   Double_t first=xmin;
 745   Double_t last=xmax;
 746   Double_t expMax=TMath::Log(last/first);
 747   for (Int_t i=0; i<nbinsX+1; ++i){
 748     (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
 749   }
 750   return binLim;
 751 }
 752
 753 //______________________________________________________________________________
 754 TVectorD* AliAnalysisTaskPIDqa::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
 755 {
 756   //
 757   // Make linear binning
 758   // the user has to delete the array afterwards!!!
 759   //
 760   if (xmax<xmin){
 761     Double_t tmp=xmin;
 762     xmin=xmax;
 763     xmax=tmp;
 764   }
 765   TVectorD *binLim=new TVectorD(nbinsX+1);
 766   Double_t first=xmin;
 767   Double_t last=xmax;
 768   Double_t binWidth=(last-first)/nbinsX;
 769   for (Int_t i=0; i<nbinsX+1; ++i){
 770     (*binLim)[i]=first+binWidth*(Double_t)i;
 771   }
 772   return binLim;
 773 }
 774
 775 //_____________________________________________________________________________
 776 TVectorD* AliAnalysisTaskPIDqa::MakeArbitraryBinning(const char* bins)
 777 {
 778   //
 779   // Make arbitrary binning, bins separated by a ','
 780   //
 781   TString limits(bins);
 782   if (limits.IsNull()){
 783     AliError("Bin Limit string is empty, cannot add the variable");
 784     return 0x0;
 785   }
 786
 787   TObjArray *arr=limits.Tokenize(",");
 788   Int_t nLimits=arr->GetEntries();
 789   if (nLimits<2){
 790     AliError("Need at leas 2 bin limits, cannot add the variable");
 791     delete arr;
 792     return 0x0;
 793   }
 794
 795   TVectorD *binLimits=new TVectorD(nLimits);
 796   for (Int_t iLim=0; iLim<nLimits; ++iLim){
 797     (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
 798   }
 799
 800   delete arr;
 801   return binLimits;
 802 }
 803